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ABSTRACT 

The prospect of detecting the first galaxies by observing their impact on the intergalactic 
medium (IGM) as they reionized it during the first billion years leads us to ask whether such 
indirect observations are capable of diagnosing which types of galaxies were most respon- 
JZLc sible for reionization. We attempt to answer this by extending the galaxy mass range of our 

first generation of large-scale radiative transfer simulations of reionization in volumes large 
enough, ^ (100/ hMpc) 3 , to make statistically meaningful predictions of observable signa- 
tures, downward from those above 10 9 Mq (high-mass, atomic-cooling halos, or "HMACHs") 
to include those between 10 8 and 10 9 Mq (low-mass, atomic -cooling halos, or "LMACHs"), 
as well. Previously, we simulated the effects of both HMACHs and LMACHs but only by re- 
ducing the box size to 35/ h Mpc, too small to apply those simulations to predict observables 
like the 21 -cm background fluctuations. Those simulations showed that LMACHS can make a 
difference, however. While LMACHs are even more abundant, photoheating suppresses their 
ability to form stars if they are located inside ionized regions of the IGM, so their contribution 
starts reionization earlier but tends to saturate over time before it ends. To predict the observ- 

■^J- ■ ables in that case while explicitly tracking the radiation from this full mass range of galaxies, 

we have had to advance both our N-body and radiative transfer methods substantially, as de- 
scribed here. With this, we perform new simulations in a box 163 Mpc on a side, and focus 
here on predictions of the 21cm background, to see if upcoming observations are capable of 
distinguishing a universe ionized primarily by HMACHS from one in which both HMACHs 
and LMACHs are responsible, and to see how these results depend upon the uncertain source 
efficiencies. 

k> , We find that 21 -cm fluctuation power spectra observed by the first generation EoR/21cm 

j_J ■ radio interferometer arrays should be able to distinguish the case of reionization by HMACHs 

alone from that by both HMACHs and LMACHs, together. Some reionization scenarios, 
e.g. one with abundant low-efficiency sources vs. one with self-regulation, yield very similar 
power spectra and rms evolution and thus can only be discriminated by their different mean 
reionization history and 21 -cm PDF distributions. We find that the skewness of the 21 -cm PDF 
distribution smoothed over LOFAR-like window shows a clear feature correlated with the rise 
of the rms due to patchiness. This is independent on the reionization scenario and thus pro- 
vides a new approach for detecting the rise of large-scale patchiness and an independent check 
on other measurements, regardless of the detailed properties of the sources. The peak epoch of 
the 21 -cm rms fluctuations depends significantly on the beam and bandwidth smoothing size 
as well as on the reionization scenario and can occur for ionized fractions as low as 30% and 
as high as 70%. Measurements of the mean photoionization rates are sensitive to the average 
density of the regions being studied and therefore could be strongly skewed in certain cases. 
Finally, the simulation volume employed has very modest effects on the results during the 
early and intermediate stages of reionization, but late-time signatures could be significantly 
affected. 

Key words: H II regions: halos — galaxies:high-redshift — intergalactic medium — 
cosmology:theory — radiative transfer — methods: numerical 
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1 INTRODUCTION 



Study of the Epoch of Reionization (EoR) has progressed in recent 
years in response to a number of new observational developments. 
The combination of the C MBR data from WMAP JKomatsu et al.l 
1201 ll ; lLarson et alj 120111) and ever deeper ground-based obser- 
vations of high-redshift QSOs, galaxies and GRBs ( Ouchi et al. 



2010 j ; |Kashikaw a et alj2o"ul;lMortlock et alj201 ll ; [Cucchiara et all 
201 lluKrug et al,l201 lb clearly suggest that the reionization process 



started early and was quite extended in time. However, observations 
of the effects of the EoR are just beginning to assemble constraints 
sufficient to diagnose the conditions which brought it about. Ongo- 
ing and upcoming observations are expected to put further, much 
more stringent constraints on the reionization history. The best 
constraints are likely to result from redshifted 21 -cm experiments 
with the low-frequency radio interferometers GMRTjjlPaciga et al. 
201 lb . LOFAF0 (e. g lHarkeretal | l2010h . MWA JLonsdale et alj 



20090 and PAPER jParsons et al]|201Ct). Additional information 
will come from Planck satellite iPlanck Collaboration et alj|201 lb 
and measurements of the near infrared background with CIBEKy 
and AKARJj, among others. 

The understanding and correct interpretation of these obser- 
vational results requires detailed modelling. Specific characteris- 
tics and features of the observable signatures provide information 
about different aspects of EoR. One of the central questions is what 
are the nature, abundances and physical properties of the ionizing 
sources. Our purpose in this work is to explore how observations of 
the EoR might diagnose the nature of the reionization sources. 

Our first generation of simulations of inhomogeneous reion- 
ization combined cosmological N-body simulations of galaxy and 
large-scale structure formation and the intergalactic density and ve- 
locity fields with detailed radiative transfer calculations of the ion- 
izing radiation from every galactic halo whose formation we re- 
solved in a volume large enough to make statistically meaningful 
predi c tions o f observable consequences of reionization (Iliev et al. 



2006 



2007 



Mellemaetai] [20 06b: Ilie v et alj|2008 a, 2007b : iDore et all 
Iliev et alj|2008bl; lHarker et alj|2009; Ichi kawa et alj|2010l : 



Fernandez et alJuOlOf ). These simulations were in comoving boxes 



of size 100//i (= 143 for h — 0.7 henceforth) Mpc on a side. Pre- 
vious simulations had been limited to much smaller volumes, too 
small to serve this purpose. Earlier simulations of smaller volumes, 
for example, underestimated the width of the time interval for the 
global tran sition of the intergalactic medium (IGM) from neutral 
to ionized (G nedin 2000l; iRicotti et alj|2002l ; lsokasian et alj|2003l ; 
ICiardi et alj 120030 , as well as the amplitude of t he kSZ fluctua- 
tions i n the temperature of th e CMB from the EOR JGnedin & Jaffd 
l200ll : ISalvaterra et alj|2005b . The characteristic size of the inter- 
galact ic H II regions during the EOR is expected to reach 10's of 
Mpc JFurlanetto et"alll2004L 12006a! : iFriedrich et alJlioTIb before 
they grow large enough to overlap. Any fluctuations introduced 
by this "patchiness" scale, therefore, require simulation volumes at 
least this large to model reionization. Moreover, since the first gen- 
eration of radio observations seeking to detect fluctuations in the 
brightness temperature of the 21cm background from the EoR have 
angular resolution of a few arcminutes, a simulation box size in 
excess of ~ 100 Mpc is necessary to characterize the power spec- 



http : //gmrt ■ ncra . tif r . res ■ in/| 

http : //www . lof ar . org/| 

http : //www . mwatele scope .org/ 

http : //physics . ucsd . edu/ -bkeat ing/CIBER. html 



trum, at wavenumbers as small as ~ 0.1 Mpc -1 where the initial 
sensitivity will peak. 

Our simulations of a comoving volume 100//i (= 143) Mpc 
on a side were limited, however, by the mass resolution of the 
N-body simulations, to the direct simulation of galactic halo 
sources more massive than ~ 2 x 10 9 AfQ. Halo sources are 
also possible at lower mass if halo gas can radiatively cool be- 
low the halo virial temperature to make star formation possi- 
ble. This includes halos above about 10 8 Mq, for which colli- 
sional excitation of H atoms can radiatively cool the primordial- 
composition halo gas since the virial temperature is above 10 4 K. 
We shall refer to these halos between about 10 8 and 10 9 Mq as 
low-mass atomic-cooling halos ("LMACHs"), to distinguish them 
from the halos above 10 9 Mq, which we shall call high-mass, 
atomic-cooling halos ("HMACHs"). LMACHs are more numer- 
ous than HMACHs at these epochs, so it might be thought that 
they would dominate reionization. However, unlike the HMACHs, 
the LMACHs are vulnerable to the negative feedback effects 
of photo-heating if they form inside a pre-exisiting H II re- 
gion of the IGM, since the pressure of the IGM would then 
prevent the intergalactic gas from collapsing gravitationally into 
their dark matter host halos (Efstathiou 1992; Shapiro et al. 199 ' 



Ouinn et aL 19961; Navarro & Steinmetz 1997; Susa & Umemur; 



ra 
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l2004llOkamotoetalJl2008l:lMesinger&Diikstrall2008b . As first 
emphasized bv lShapiro et alj ( 11994b , this limits their contribution 
to reionization. Halos of even smaller mass than LMACHs would 
be even more vulnerable to the negative feedback effects of pho- 
toionization heating, since they would, in addition to being pre- 
vented from capturing intergalactic gas, photoevaporate whatever 
interstellar gas they ha d already a ccumulate d before they wer e en- 
gulfed by reionization dShapiro e t al. 2004; Iliev et al. 20051). Be- 
fore that reionization, however, these minihalos with masses below 
about 10 8 Mq, nevertheless could have formed stars if enough H2 
molecules were present in them. Since their virial temperatures are 
below 10 4 K, their gas is too cold for collisional excitation and ra- 
diative cooling by H atoms, but cooling through the collisionally 
exciting rotational-vibrational levels of H2 is possible. However, 
these molecules were easily dissociated by the rising UV back- 
ground of starlight at energies below the ionization threshold of H 
atoms, an inevitable by-product of the same stars that contributed 
to reionization. This tends to limi t the contribution of minihalos 
to rei onization early in the EoR JHaiman et al] l2000l : lAhn et al] 
2009). Previous estimates of the minihalo contribution, as a result, 
assume that this contribution is smaller than that of the LMACHs 
and HMACHS, so for now, we shall neglect it, although patchiness 
in the UV ba ckground may m ake them important than one might 
naively think JAhn et alj2009b . 

To investigate the impact of the LMACHs on reionization 
and its observable properties by direct radiative transfer simula- 
tion of reionization, our first generation of simulations boosted 
the halo mass resolution in order to resolve all the halos of mass 
10 8 Mq and above, but sacrificed volume by simula ting in a box 
of size 37 /h (= 53 Mpc) on a side (I liev et alj|2007ab . These sim- 
ulations demonstrated explicitly that reionization in the presence 
of the LMACHs and their suppression if they formed inside pre- 
existing H II regions during the EoR was "self-regulated". The 
more LMACHs that formed, the more volume and mass of the IGM 
was ionized, but as this ionized fraction grew, so did the fraction of 
the total LMACH halo population that formed inside the ionized 
regions and was suppressed as sources of reionization. This meant 
that, although the LMACHs dominated the early phase of reion- 
ization, their contribution to reionization eventually saturated, and 
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reionization was finished by the HMACHs, whose abundance rose 
exponentially over time. An observational consequence of impor- 
tance that is a possible signature of this process was an early onset 
but late finish for the EoR, as required to explain the high values 
of electron scattering optical depth reported from WMAP observa- 
tions of the large-angle fluctuations in the polarization of the CMB. 
To make further predictions of observable consequences, however, 
we need to enlarge the volume of these simulations to greater than 
100/h Mpc on a side, without losing this enhanced mass resolution 
necessary to resolve the LMACHs. That is the purpose of the new 
developments we report in this paper. 

To accomplish this goal, it was first necessary to improve and 
advance both our N-body and radiative transfer methods in order to 
be able to simulate halo formation with much higher mass resolu- 
tion and to transfer the ionizing radiation from a much larger num- 
ber of sources. Toward this end, both codes had to become mas- 
sively parallel, running on thousands of computing cores, as well 
as more efficient. These nume rical developments are discussed in 
more detail in §[2]below and in llliev et al.l J2008c). 

While the work described here follows naturally from our own 
previous work as described above, it also differs substantially from 
other work in the literature to-date involving large-scale radiative 
transfer simulations of reionization. A full account of that literature 
is well-beyond the scope of this paper, but we will mention a few 
points to distinguish the current work. Our N-body simulations re- 
solve all galactic halo sources of 10 8 Mq and above, in a comoving 
box as large as lU/h = 163 Mpc on a side. We post-process the 
density field of the IGM and the galactic halo source populations 
derived from these N-body simulations by performing a detailed, 
ray-tracing calculation on a grid of 256 3 cells. The simulations de- 
scribed in McOuinn et al. (2007) were based upon post-processing 
N-body simulations with halo mass resolution above 10 9 Mq, in 
a box 65.6/h Mpc on a side, a volum e which is five times s maller 
than ours, also on a grid of 256 3 cells. iMcOuinn et alj ( 120070 were, 
thus, unable to treat explicitly the halo mass range below 10 9 Mq, 
which is subject to suppression by the negative feedback effects 
of reionization, but they did include a semi-analytical, "subgrid" 
approximation for the contribution from the unresolved, smaller - 
mass halos, including some feedback effects. iTrac & Cerj (2007) 
considered a 50/ h Mpc box, hence, an order of magnitude smaller 
volume than ours, and 180 3 cells for radiative transfer, but their 
halo mass resolution was similar to ours. They did not consider 
the feedback effects of reionization on t he small-mass galactic halo 
sources as we do here. lShinet alJ (2008) subsequently applied this 
method to simulate a volume 100/ft Mpc on a side (2/3 of our 
volume), with similar halo mass resolution, on a radiative transfer 
grid of 360 3 cells, to study the structure of the patchy ionization 
field during the EOR, again without considering feedback and do- 
ing only a single simulation, without considering any variation of 
th e (highly uncertain) rei oniza tion parameters. Re cently, the codes 
of IMcOuinn et al.l <2007h and Trac & Cenl 12007) were compared 
with each other in IZahn et alJ q201 lh . by comparing results for 
the ionization fields and 21 -cm brightness temperature fluctuation 
statistics, for reionizaton simulations advanced part-way through 
the epoch of reionization, to the 72% ionized point, based upon 
post-processing a previously-simulated input density field in a box 
100/h Mpc on a side, smoothed to a radiative transfer grid with 
256 3 cells, with halo mass resolution of 10 8 solar masses, without 
any feedback or back-reaction on the halo sources. They also con- 
sidered a single reionization scena rio, with no parameter variation. 
Finally, lAubert & Tevssiej ( 1201 Of) simulated radiative transfer in 
several different boxes, as large as 100/h Mpc, by post-processing 



a density field and galaxy population derived from separate sim- 
ulations that combined N-body dynamics and hydrodynamics on 
a 1024 3 grid, but with minimum resolved galaxy masses in that 
case as large as 8 x 10 9 Mq. No feedback from reionization on 
galactic sources was considered. While there are other distinctions 
of detail both between these other simulations and ours and of one 
from another, we have listed these above to make it clear that our 
current paper will describe simulations and results which are new, 
both because they are based upon different methodology, applied in 
greater depth than previously to predict observables like the 21cm 
background from the EOR, and because they are on an unprece- 
dented scale. 

The rest of paper is organized as follows. In § [2] we present 
our codes, numerical methods and simulations. In §[3] we discuss 
our results on the formation of early cosmic structures. In §[4] we 
present our results on the basic reionization features, reionization 
history, integrated electron-scattering optical depth and geometry 
of ionized patches. The observational signatures derived from our 
simulations are presented and discussed in §[5] Our conclusions are 
summarized in § [6] Finally, in Appendix [A] we present the (physi- 
cally less realistic) cases of reionization by rare, massive sources, 
while in Appendix [B] we discuss the more technical point of the 
dependence of our results on the Jeans suppression threshold for 
low-mass sources. 



2 SIMULATIONS 



Our ba s ic methodology has bee n pr eviously de s cribed in llliev et al.1 
d2006h : iMellema et all J2006b1) and llliev et all J2007ah . Due to the 
much larger scale of our current simulations compared to our pre- 
vious ones, both our structure formation and our radiative transfer 
code had to be significantly developed and re-designed, in partic- 
ular to allow their massive paralellization on distributed-memory 
machines. In this section we present our new set of simulations, 
along with a summary of our methods and parallel code scaling to 
large number of computing cores. 



2.1 N-body simulations 

We start by performing very high resolution N-body simulations 
of the formation of high-redshift structures. We use the CubeP 3 M 
N-body codqj which evolved from the particle-mesh (PM) code 
PMFAST rMerzetalJ (2005). In CubeP 3 M several important new 
features were introduced in comparison with these previous codes. 
The first one is the addition of a short-range direct particle- 
particle force, making it a P 3 M (particle-particle-particle-mesh) 
code. This significantly improves its spatial resolution and accu- 
racy at small scales compared to PM codes. A second important 
new development is that CubeP 3 M is now a massively parallel 
code which can run efficiently on either distributed- or shared- 
memory machines. This is achieved through cubical equal-volume 
doma in decomposition and a hybrid MPI and OpenMP approach 
(see dlliev et al] l2008ch for more details). CubeP 3 M scales well 
(with 'weak' scaling, whereby the execution time rises proportion- 
ally to the problem size) up to thousands of processors, as shown 
in Fig. Q] (left) and has to date been run on up to 21,976 comput- 
ing cores, following up to 5488 3 particles dlliev et alfcOlOO . These 



6 |http : //www . cit a . utoronto . ca/mediawiki/index . php/CubePM 
for description of the code see also dlliev et al. 2008c). 
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Table 1. N-body simulation parameters. Background cosmology is based on the WMAP 5-year results. 



boxsize 



Npart mesh spatial resolution ^particle. 



M ha i . 



37 h^Mpc 1024 3 2048 3 1.81 kpc/h 5.05 x 10 6 M 1.01xl0 8 M Q 

64 /i _1 Mpc 1728 3 3456 3 1.85 kpc/h 5.44 x 10 6 M Q 1.09 x 10 8 M Q 

74 h-^-Mpc 2048 3 4096 3 1.81 kpc/h 5.05 x 10 6 M Q 1.01xl0 8 M Q 

114 h- 1 Mpc 3072 3 6144 3 1.86 kpc/h 5.47 x 10 6 M Q 1.09xl0 8 M Q 



scaling tests were run on the Texas Advanced Computing Center 
(TACC) computer, Lonestar, and on the currently-available portion 
of the European Petascale computer under development in France, 
Curie, at CEA, part of the Partnership for Advanced Computing 
in Europe (PRACE). Our results show almost perfect scaling of 
CubeP 3 M, within 3% of the ideal one (dashed line), for up to 2,048 
cores. 

We performed a series of N-body cosmic structure forma- 
tion simulations (summarized in Table QJ with increasing simu- 
lation box size, from 37/i _1 Mpc (53 Mpc) up to 114ft _1 Mpc 
(163 Mpc), but all with a fixed spatial and mass resolution. These 
N-body simulations were run on a range of core number from 256 
(1024 3 particles) up to 2048 cores (3072 3 particles). The force 
smoothing length is fixed to 1/20 of the mean inter-particle spac- 
ing, or 1.8 h^ 1 kpc. The largest of these simulations follows a vol- 
ume which is 50% larger than the largest structure formation simu- 
lation performed previously at similar resolution. 

The N-body simulations required between 4,100 (for 1024 3 
particles) and 159,000 (for 3072 3 particles) computing hours (com- 
puting cores x wall-clock hours) on the TACC computer Ranger 
(SunBlade x6420 with AMD x86 64 Opteron Quad Core, 2.3 GHz, 
9.2 GFlops per core Barcelona processors and Infiniband network- 
ing). The particle mass is 5 x 1O 6 M0, which guarantees that all 
atomically-cooling halos (M > 10 8 Mq) are resolved with at 
least 20 particles. We use a spherical overdensity halo finder with 
overdensity parameter fixed to 178 of mean density (Note that 
pM 7^ Pcrit except for very high redshift. The series of N-body 
simulations with an increasing size allowed us to test the conver- 
gence of our results with computational box size. 

The background cosmology is based on WMAP 5-year data 
combined with constraints from baryonic acoustic oscillations 
and high-redshift supernovae (Qm — 0.27, £7a = 0.73, h = 
0.7, tt b - 0.044,0-8 = 0.8, n - 0.96). The linear power spec- 
trum of density flu ctuations was calculated with the code CAMB 
JLewis et alj uOOOl) . Initial conditions were generated using the 
Zel'dovich approximation at sufficien tly high redshift ( zj = 300) 
to ensure against numerical artifacts JCrocce et alj|2006r) . 



2.2 Radiative transfer simulations 

The radiative transfer simulations are perfo rmed with our code C 2 - 
Ray (Conservative Causal Ray-Tracing) JMellema et alj 12006a). 
The method is explicitly photon-conserving in both space and time 
for individual sources and approximately (to a good approxima- 
tion) photon-conserving for multiple sources, which ensures cor- 
rect tracking of ionization fronts without loss of accuracy, indepen- 
dent of the spatial and time resolution, with corresponding great 
gains in efficiency. The code has been tested in detail ag ainst a num- 
ber of exact analytical solutions (Mellema et al. 2006a), as well as 



in direct comparison with a number of other independent radia- 
tive transfer methods on a standardiz ed set of benchmark problems 
dlliev & et a l. 2006; Ili ev et al.l2009l) . The ionizing radiation is ray- 
traced from every source to every grid cell using the short char- 
acteristics method, whereby the neutral column density between 
the source and a given cell is given by interpolation of the column 
densities of the previous cells which lie closer to the source, in 
addition to the neutral column density through the cell itself. The 
contribution of each source to the local photoionization rate of a 
given cell is first calculated independently, after which all contri- 
butions are added together and a nonequilibrium chemistry solver 
is used to calculate the resulting ionization state. Ordinarily, multi- 
ple sources contribute to the local photoionization rate of each cell. 
Changes in the rate modify the neutral fraction and thus the neutral 
column density, which in turn changes the photoionization rates 
themselves (since either more or less radiation reaches the cell). 
An iteration procedure is thus called for in order to converge to 
the correct, self-consistent solutio n. While our basic m ethodology 
remains essentially as described in lMellema et alj d2006af) . our C 2 - 
Ray code has been thoroughly re-written in Fortran 90, made more 
flexible and modular and parallelized for distributed-memory ma- 
chines. In terms of parallelization strategy, due to the causal nature 
of the ray-tracing procedure (i.e. the state of each cell can be cal- 
culated only after all previous cells, closer to the source are done) 
it is not possible to employ domain decomposition (except for a 
limited one, into octants, see below), although other approaches 
exist which seek ways to ov ercome this limitation dNakamoto et alj 
1200 lc iRiikhorst et alj |2006). Instead, the main code loop over the 
sources of ionizing radiation is done in massively parallel fash- 
ion. Each MPI node has a copy of the density field and receives 
a number of sources whose radiation is to be traced through the 
grid. For the large-scale cosmological reionization problem there 
are typically hundreds of thousands to millions of sources, thus 
our code scales well up to tens of thousands of cores at least (see 
next section). For problems with (relatively) low number of ioniz- 
ing sources such parallelization strategy would be inefficient, but 
such problems are not sufficiently computationally-intensive to re- 
quire such massive parallelization and could, instead, be solved on 
a smaller number of nodes, or even in serial. A similar situation oc- 
curs for the initial steps of the simulations presented below, when 
the cosmological structure formation is not yet much advanced, 
thus only a few to few tens of halos form. However, their num- 
ber increases exponentially over time, quickly reaching thousands, 
and then tens and hundreds of thousands. We therefore start our 
simulations on a small number of cores (typically 32), raising to 
thousands of cores as more sources form. 

As mentioned above, a limited domain decomposition onto 
octants is possible for our method, since those are independent of 
each other within the short-characteristic ray-tracing framework. 
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Figure 1. (a)(left) Scaling of the CubeP 3 M code. Plotted are the code speedup vs. the number of computational cores used. Both quantities are normalized to 
the smallest run in each case. (b)(right) Scaling of the C 2 -Ray code. Plotted are the code speedup vs. the number of computational cores used, again normalized 
to the smallest of the three runs compared. Dashed line indicates the ideal weak scaling for each case. 

Table 2. Reionization simulation parameters and global reionization history results. All runs use background cosmology based on the WMAP 5-year results. 
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run 


boxsize 
[cMpc] 
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mesh 
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[Mq] 


min unsupp. 
halo [Mq] 
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Z W% 
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3 f 1 is related to g 7 by Eqn. with At = 11.53 Myrs. 



We use this to (optionally) improve the memory efficiency of the 
code by doing the grid octants in parallel within each MPI node 
using OpenMP multi-threading. This way each MPI node needs 
only one copy of the grid, which is shared amongst the cores within 
the node. 

The radiative transfer problem size scales proportionally to 
both the grid size and the number of sources. Results, shown in 
Figure [T] (right) demonstrate almost perfect scaling, within ~ 10% 
from the ideal one, for up to 8,192 cores. 

The N-body simulations discussed above provide us with the 



spatial distribution of cosmological structures and their evolution 
in time. We then use this information as input to a full 3D ra- 
diative transfer simulations of the reionization history, as follows. 
We saved series of time-slices, both particle lists and halo cata- 
logues from redshift 50 down to 6, uniformly spaced in time, every 
At = 11.53 Myr, a total of 76 slices. Based on the particle distri- 
bution at each redshift we used SPH-style smoothing scheme using 
the nearest neighbours (to be described in detail in a companion 
paper, in prep.) to produce regular-grid density and bulk velocity 
fields at the radiative transfer resolution of 256 3 cells. 
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All identified halos are potential sources of ionizing radiation, 
with a photon production rate per Myr, N 1 , proportional to their 
mass, M: 



N-, 



f-yMtlb 

At (iQarrip 



(1) 



where m p is the proton mass and / 7 = f CS cf*N- k is an ioniz- 
ing photon production efficiency parameter which includes the ef- 
ficiency of converting gas into stars, /*, the ionizing photon es- 
cape fraction from the halo into the IGM, / csc and the number of 
ionizing photons produced per stellar atom, N* . The latter param- 
eter depends on the assumed IMF for the stellar population and 
varies between 4,000 and ~ 100, 000. Halos were assigned differ- 
ent luminosities according to whether their mass was above ("large 
sources") or below ("small sources") lO 9 Af0 (but above 1O 8 M0, 
the minimum resolved halo mass). Low-mass sources are assumed 
to be suppressed within ionized regions (for ionization fraction 
higher than 10%) , through Jeans-mass filtering, as discussed in 
llliev et al.l J2007ah . 

We note that while previously we used the factor, / 7 , to char- 
acterize the source efficiencies, here we define a slightly different 
factor, g~f, that is given by 



fh 



=/r(^) 



(2) 



where At is the time between two snapshots from the N-body sim- 
ulation. The new factor g 1 has the advantage that it is independent 
of the length of the time interval between the density slices, and 
as such it allows a direct comparison between runs with different 
At. For reader's convenience we listed the values of both parame- 
ters in Table [2] We also note that the specific numerical values of 
the efficiency parameters are strongly dependent on the background 
cosmology adopted and the minimum source halo mass. Therefore, 
parameter values for simulations based on different underlying cos- 
mology and resolution should not be compared directly, but would 
require a cosmology and resolution-dependent conversion coeffi- 
cients to achieve the same reionization history. 

Our full simulation notation reads Lbox_gI_J(S)(K) (the 
bracketed quantities are listed only when needed), where ' Lbox' 
is the simulation box size in Mpc, 'J' and 'J' are the values of the 
ffgamma factor for HMACHs and LMACHs, respectively, the sym- 
bol 'S' means that the small sources are suppressed within already- 
ionized regions and 'K' indicated the ionized fraction threshold for 
a given radiative transfer cell above which this suppression occurs 
for halos residing in that cell, which is 0. 1 if not listed explicitly 
and raised to 0.9 or 0.5 for cases K — 9 and K = 5, respectively 
(see below for details). For example, 53Mpc_g8.7_130S indicates 
that large sources have an efficiency g 1 — 8.7, while small sources 
have an efficiency <? 7 = 130 and are suppressed in ionized regions. 

We have performed series of radiative transfer simulations 
with varying underlying assumptions about the source efficiencies 
and the suppression conditions imposed on the low-mass sources, 
as summarized in Table[2] For our radiative transfer simulations we 
use the data from the largest N-body box, 114/h = 163 Mpc, and 
the smallest one of these, 37//i = 53 Mpc. The former volume 
is sufficiently large to faithfully represent the reionization observ- 
ables, while the latter one affords much faster and computationally 
cheaper simulations, which allows us to explore a wider parameter 
space. These two very different computational volumes also allow 
us to investigate resolution effects and evaluate which features of 
reionization and observable signatures are sensitive to the box size 
and which are less so. We label all runs by a short label (listed in 



the first column of Table[2} for more compact notation. Large-box 
runs are labelled L1-L3, while small-box ones are labelled S1-S9. 

Our fiducial runs, to which all the others will be compared are 
163Mpc_g8.7_130S (LI) and the companion small-box one with 
same source efficiencies, 53Mpc_g8.7_130S (SI). These parame- 
ters yield a relatively early overlap and high electron scattering op- 
tical depth. The second set of simulations, 163Mpc_gl.7_8.7S (L2) 
and 53Mpc_gl.7_8.7S (S2) are in the opposite limit, which assumes 
considerably lower efficiencies for both types of sources and as a 
consequence serves as a model for a late-overlap, extended, more 
photon-poor reionization scenario. These two cases are designed to 
roughly bracket the range of observationally-allowed reionization 
scenarios. 

Our third large-box simulation, 163Mpc_g21.7_0 (L3) is 
equivalent to our prev ious large-box si mulations without self- 
regulation presented in jlliev et al.ll2008aO . except for the updated 
background cosmology and the source photon production efficies, 
adjusted here to yield the same overlap epoch as our fiducial case 
LI. Therefore, LI and L3 share the underlying density structures 
and sources (apart from the different minimum mass cutoffs) and 
hence a head-to-head comparison yields the effects of the presence 
of low-mass sources and Jeans-mass filtering. 

The rest of our cases, S3 to S9, test various aspects of the 
reionization source modelling. Simulation S3 is an extreme case 
which has the same source efficiencies as our fiducial case LI, 
but assumes no suppression occurs. Naturally, this results in a 
very early reionization and very high integrated optical depth, 
t cs = 0.131, which is well outside the WMAP5 1 — a range of 
t os = 0.084 ± 0.016. Simulation 53Mpc_g0.4.5.3 (S4) is again 
without suppression, but here we tuned down the efficiencies of 
both types of sources so as to achieve approximately the same over- 
lap epoch as in our fiducial case LI. In simulation 53Mpc_gl0.4_0 
(S5) we assume that there are no low-mass, M < 1O 9 M0 sources 
at all and we adjust the photon efficiency of the remaining sources 
to again reach overlap at roughly the same epoch as in the fiducial 
case. 

Additionally, we consider two scenarios which have exactly 
the same time-dependent ionizing photon emissivity (and therefore 
almost identical reionization history) as our fiducial case SI, but 
with higher minimum source mass of 1O 9 M0 (53Mpc_uvS_le9; 
S8) and 10 10 M Q (53Mpc_uvS_lel0; S9). The fixed ionizing pho- 
ton emissivity results in unphysically high early source luminosi- 
ties and we consider them primarily in order to illustrate the ef- 
fect of re-distributing the full luminosity of all sources over the 
massive on es only, similar to the models adop ted in some re- 
cent work JThomas et alj|2009l : iBaek et al. 2009). These simula- 
tions and some illustrative results from them are discussed in Ap- 
pendix [A] 

Finally, in order to evaluate the rubistness of our 
source suppression model, we consider two more scenarios, 
53Mpc.g8.7.130S9 (S6) and 53Mpc.g8.7.130S5 (ST), whereby 
we raise the ionization threshold for low-mass source suppression 
to ^threshold = 0.9 and 0.5, respectively, from our fiducial 
threshold of ^threshold — 0-1- Since this is a more technical study 
we present its results separately, in Appendix iBl 

The radiative transfer simulations presented in this work typi- 
cally required ~ 0.5 — 1 million computing hours (163 Mpc boxes) 
and ~ 10 — 30 thousand computing hours (53 Mpc boxes), depend- 
ing on the specific set of source parameters we adopted. 
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Figure 2. (left) Slice of the Cosmic Web at redshift z = 6 from our CubeP M simulation with 3072 particles (29 billion) on a 6144 fine grid in a comoving 
volume of 163 Mpc on a side. Shown are the dark matter density (blue) and halos (in actual size; yellow). Image resolution is 6144 X 6144, the slice is 1/h 
Mpc thick, (right) Zoomed-up region (25.76x25.76 Mpc) of the same image. 
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Figure 3. (left) Halo abundances for v = 6 c D + /a(0, M)=l (M„; 
Collapsed fraction of HMACHs (top) and LMACHs (bottom) halos. 



black, solid), 2 (red, dotted), 3 (green, short-dashed), and 5 (blue, long-dashed), (right) 



3 RESULTS: EARLY STRUCTURE FORMATION 

In Figure [2] (left: full box, right: zoomed sub-volume) we show 
a slice of the density field and halos at redshift z — 6 from our 
114 h' 1 Mpc (163 Mpc), 3072 3 -particle N-body simulation. The 
structure formation is already well-advanced and strongly nonlin- 
ear at sub-Mpc scales. The very first resolved (Afhaio > 10 8 Mq) 
halo in this volume forms at z — 31, while the first insupppressible 
halo (A/haio > 10 9 Mq) forms at z = 21. By z — 6 there are over 
20.5 million collapsed halos, of which ~ 18.7 million low-mass 
(Mhaio < 1O 9 M ) halos and ~ 2 million high-mass (Afhalo > 



10 Mq) halos. The halos are strongly clustered at all times, more 
so going to higher redshifts, when they are ever rarer. The halo 
abundances are usually described in terms of v — S c D+/a(0, M), 
where 5 C is the linear density contrast corresponding to the moment 
of collapse of a top-hat density perturbation, <r(0, M) is the present 
variance of the density fluctuations corresponding to the mass scale 
M, and D+ is the growth factor of the density fluctuations. The 
halo masses corresponding to v — 1 (most common, M*, halos), 
v — 2, 3, (rare halos) and v — 5 (extremely rare halos) are shown 
in Figure[3](left). Clearly, before ~ 6, within the redshift range of 
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Figure 4. Simulated halo mass function at high-z derived from 1 14/h Mpc box (solid, red), 74/h Mpc box (dot-dashed, green), 37/h Mpc box (long-dashed, 
blue) at (left to right) z = 12, 9 and 6. Also shown are the Press-Schechter (dotted, black) and Sheth-Tormen (short-dashed, black) analytical mass functions. 
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Figure 5. The power spectra of the density fields, A p = (fc 3 P(fc)/2-7r 2 ) 1 / 2 , for the 114/fe = 163 Mpc box run (red, solid) and the 37/h = 53 Mpc box 
run (blue, dashed) at redshifts z = 12 (left) and 2 = 9 (right). 



interest here, there are no 1 — a halos at all and all halos are rare. 
Atomically-cooling halos are 2 — 3 — a at the low end of the redshift 
interval and as rare as 5 — a at early times. The evolution of the col- 
lapsed fractions in high-mass and low-mass halos is shown in Fig[5] 
(right). The collapsed fractions start very low and rise exponentially 
at early times when the halos are very rare. The collapsed fraction 
in low-mass halos reaches 1CP 3 and 1% at z = 14 and z = 9.8, 
respectively. After that point it starts to level off as low-mass halos 
start to become less rare (2 — a or less) and their collapsed fraction 
reaches 3.4% by z — 6. The collapsed fraction in high-mass ha- 
los rises steeply all the way to z = 6, eventually reaching 4.25%. 
There is only a modest departure from the exponential growth, re- 
flecting the fact that they remain quite rare throughout this period. 
The simulation volume has essentially no effect on the derived col- 
lapsed fractions, indicating numerical convergence on that quantity 
at fixed mass resolution. The only exception to this is at very early 
times (z > 26 for the low-mass halos and z > 17 for the high- 



mass ones) the corresponding halo populations are so rare (~5-<r 
in each case) that Poisson noise (i.e. cosmic variance, due to the 
smaller volume) affects the results. For example, the first resolved 
halos form at Z — 31 in the 114 /i -1 Mpc box, but only at z — 26 
in the 37/i _1 Mpc box. As soon as there is sufficient statistics for 
any given volume the collapsed fractions converge. 

The halo mass functions derived from our simulations are 
shown in Figure [4] for a ran ee of redsifts, z = 12 — 6, along 
with the Pre ss-Schechter (PS; Press & Schechtetll 19741) and Sheth- 
Tormen (ST; Sheth & Tormen 2002) analytical mass functions. The 
halo abundancies at all redshifts fall between those two analytical 
predictions. PS always under-predicts the abundances of massive, 
rare halos, while ST over-predicts them. With time ST becomes a 
better match to the numerical results. This broadl y agrees with pre- 
vious results on t he high-redshift m ass functions Jlliev et alJI 20061 ; 
iReed et alj2"007| ; |Lukic et alj|2007t) . 

In Fig. [5] we plot the total matter density field power spectra 
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Figure 6. The halo bias, 6 hh = A h h/A p , for the 114//i = 163 Mpc box run (red, solid) at redshifts z = 12.04 (left) and z = 9.03 (right). Lines are for 
halos binned by decades of mass (bottom to top curve) 1O 8 M < A/ halo < 1O 9 M , 1O 9 M < A/ halo < 1O 1O M , 10 10 Af Q < A/ halo < lO n M , 
and lO u Af < Afhaio < 1O 12 M . 
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Figure 7. The effect of self-regulation on the reionization history and integrated electron-scattering optical depth: (left) Mass-weighted reionization histories 
(bottom) and the ratio of the mean mass-weighted and volume ionized fractions, x m /x v (top) for our fiducial self-regulated case, LI (blue, solid) and the 
corresponding non-selfregulated case with same overlap epoch, L3 (black, long-dashed). The computational box size is 163 Mpc in both cases. Inset shows 
the same reionization histories in linear scale, (right) The corresponding electron scattering optical depth, t cb (z) integrated from redshift to redshift z for 
the same two simulations. Horizontal lines indicate the mean and 1-cr band derived from the WMAP 5-year data, while the dotted line shows the value of r es 
for a fully-ionized universe. 



A p = (k 3 P(k)/2ir 2 ) 1 ' 2 at two representative redshifts for our 
largest (114 Mpc/h) and smallest (37 Mpc/h) boxes. Power spec- 
tra were calculated by interpolating the N-body particles using a 
cloud-in-cell scheme onto the fine grid of the CubeP 3 M code, with 
6144 3 and 2048 3 cells, respectively. There is a close agreement 
between the two cases, apart from the expected variance at scales 



close to the box size. This shows that there is no missing density 
fluctuation power in the small box, except for the scales at or above 
the box size. 

Finally, in Figure[6]we show the halo-halo bias, calculated as 
the ratio of the halo autocorrelation power specrum divided by the 
density field one, i.e. bhh = Ahh/A p . This measures the cluster- 
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Figure 8. (bottom) Number of ionizing photons emitted by all active 
sources in the computational volume per timestep and (top) cumulative 
number of photons per total gas atom released into the IGM. Notation 
is the same as in Fig. [JJ The vertical line marks the overlap redshift 
(z(x m = 0.99)) for each case. 



ing of the dark matter halos with respect to the underlying mat- 
ter density field at redshifts z = 12 and 9, roughly corresponding 
to early and advanced stages of reionization. Because of the rela- 
tive rarity of all halos studied here, their clustring is quite strong, 
particularly at the smallest resolved scales (k ~ 100 h/Mpc), 
where it is of order 100 for the lowest mass halos and is as high 
as ~ 10 4 for the most massive halos. At these small scales the 
bias is strongly nonlinear. At large scales the bias factors asymp- 
tote to a (m ass-dependent) constant - the large-scale linear bias 
I IMo & Whi te 1996) The transition between the large-scale linear 
bias and the nonlinear one occurs around k ~ 1 for the lowest- 
mass halos at z — 9, rising to k ~ 0.1 — 0.5 for the larger halos 
and/or higher redshifts. The largest, rarest halos (M > 1O 1O M0 at 
z — VIM > 10 11 M Q at z — 9) show a roughly linear log-log 
b(k) relation and never asymptote to the linear bias value within 
the k-range covered in our simulation. 



4 RESULTS: BASIC FEATURES OF THE 
SELF-REGULATED REIONIZATION 

4.1 The effects of self- regulation 

We start by comparing our fiducial simulation LI against the sim- 
ulation L3, equivalent to our previous l arge-box simulati ons with- 
out LMACHs and their self-regulation dlliev et alfeoOSal) . The re- 
sulting reionization histories and integrated electron-scattering op- 
tical depths are shown in Figure [7] As expected, the self- regulated 
model yields a much more gradual and extended reionization his- 
tory, which starts with the formation of the first atomically-cooling 
halo sources at z ~ 31 vs. the much later start, at z ~ 20 in L3, due 
to the much larger, rarer sources in the latter case, which accord- 
ingly form later. The exponential rise in the numbers of the high- 
mass sources yields a steep, power-law like reionization history 
(reasonably well-fit by lg(z m ) « -1.226,2+10.41 for a; m > 0.03, 



rising somewhat steeper than this earlier on) when it is driven solely 
by those sources, while in the self-regulated case the steep ini- 
tial rise becomes much more gradual when the self-regulation first 
kicks in, around z ~ 16, when the efficient low-mass sources are 
massively suppressed and thereby gradually give way to the less ef- 
ficient high-mass ones which come to dominate at the latter stages 
of reionization. However, we note that even with self-regulation the 
reionization history remains monotonic and no plateaus, let alone 
double reionization, ever occur. The mass-weighted over volume 
ionized fraction (upper panel) is always lower in the self-regulated 
case, indicating that reionization has less pronounced inside-out 
character, i.e. ionized regions are less correlated with the highest 
density peaks in this case since reionization is driven by wider 
range of sources, including low-mass, less biased ones. The in- 
tegrated electron-scattering optical depth (Figure [7] right) is sig- 
nificantly boosted by the presence of low-mass sources, by about 
0.01 overall, most of it due to the early stages of reionization. For 
the particular source efficiencies we have chosen here both optical 
depths fall within the 1-ct interval given by the WMAP 5-year data, 
albeit the value for the self-regulated fiducial case is very close to 
the central value, while the r cs for the L3 case is at the low 1-ct 
limit. 

These reionization histories are a direct consequence of the 
overall number of ionizing photons emitted by all active sources, 
shown in Figure[8] In the case L3 where no source suppression oc- 
curs the number of photons emitted per timestep simply rises pro- 
portionally to the halo collapsed fraction, roughly exponentially. 
In contrast, in the fiducial self-regulated case LI the initial expo- 
nential rise is halted around redshift z ~ 16 and rises very slowly 
(and moderately non-monotonically) until z ~ 11, at which point 
sufficient number of high-mass, non-suppressible sources form to 
allow them to take over the evolution, while the low-mass sources 
becom e highly suppressed . Therefore, similarly to our earlier re- 
sults in llliev et al.l J2007af) the late phase of reionization and over- 
lap epoch, z ov , are dominated by HMACHs, while the LMACHs 
dominate the early phase of reionization and provide a significant 
boost to the electron-scattering optical depth, r cs . Ultimately, by 
overlap in both simulations LI and L3 there are 1.2-1.6 ionizing 
photon per atom emitted, slightly more in the self-regulated case 
due to its more extended reionization history which yields more 
recombinations per atom. 

4.2 The effects of source efficiencies and box size 

The reionization histories derived from our suite of simulations are 
shown in Fig. [9] The reionization history is monotonic in all cases, 
although due to the self-regulation the slope of the curves can vary 
significantly and in particular can become almost horizontal for 
short periods of time when the Jeans mass filtering compensates 
for the rise in source numbers. The exact redshifts at which certain 
reionization milestones, 10%, 50% and 90% by mass, are reached 
are listed in Table[2] We also list there the epochs when final over- 
lap, which we define as the time when x m — 0.99, i.e. at least 99% 
of the mass is ionized, is reached in each case. 

Our large- volume, self-regulated simulations LI and L2 
(Fig.|9l left panels) have reionization histories which are very sim- 
ilar to each other, but offset by Az ~ 2. Overlap is reached 
at z — 8.3 (6.7) in LI (L2), corresponding to early (extended) 
reionization scenarios. The integrated electron scattering optical 
depth for LI is r cs = 0.080 for the early reionization case, well 
within the current WMAP5 1-ct constraints. The corresponding 
value for the low-efficiency, extended reionization scenario L2 is 
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Figure 9. (bottom panels) Mass-weighted reionization histories for our fiducial self-regulated cases (left) and with varying assumptions about the ionizing 
sources and their suppression (right), (top panels) Ratio of the mean mass-weighted and volume ionized fractions, x m /x v . All cases are labelled by color and 
line-type, as follows: (left) LI (blue, solid), SI (green, short-dashed), L2 (magenta, long-dashed), and S2 (red, dotted), (right) S4 (cyan, long dashed), S5 (light 
red, dot-short dashed), and S3 (light green, dot-long dashed). For ease of comparison we show the fiducial case S 1 on both plots. 



0.058, which is outside the 1-cr range, but still within 2-cr. On 
the other hand, the simulation volume (163 vs. 53 Mpc) has lit- 
tle effect on the global reionization histories. This is in agree- 
ment with the results in Iliev et al. (2006), which were derived by 
sub-dividing a 100 /i _1 Mpc volume into smaller ones, which indi- 
cated that ~ 20 — 30 ft -1 Mpc box is sufficient to reliably derive 
the global mean reionization history. Most variations between the 
corresponding large and small box simulations result solely from 
the different random realizations in the two cases. At early times 
(x m < 0.01) there are also departures due to cosmic variance - un- 
like the larger, 163 Mpc, volume the more limited 53 Mpc one does 
not contain any sources at z > 25 as those are statistically too rare 
to occur. Even when the very first halos appear in the 53 Mpc vol- 
ume, they are initially so few that they are subject to very high shot 
noise fluctuations. Once there are statistically-significant numbers 
of sources in each size box the reionization histories converge and 
any fluctuations thereafter are simply due to the different random 
realizations. There is also some effect from the higher resolution of 
the small-box simulations, due to the better-resolved density field 
in those cases, which yields slightly increased recombinations. This 
effect is rather minor here however, because the relatively small dif- 
ference in resolution results in only a marginal increase of the the 
recombination rates. 

In Fig.|9](top left panel) we show the ratio of the ionized frac- 
tion by mass, x m and by volume, x v , which is equal to the aver- 
age d ensity of the ionized regions in units of the mean (Ili ev et al] 
2006). These ratios start at about 2 and remain above unity at all 
times, indicating that the reionization proceeds in an inside-out 
manner, with the high density peaks (where the first sources pref- 
erentially form) being ionized on average earlier than the mean and 
low-density ones. On average the ionized regions are denser in the 
low-efficiency cases. This behaviour could be expected based on 
the typically smaller H II region sizes in those cases. They there- 
fore stay in the immediate vicinity of the density peaks and do not 



propagate as much into the voids. The higher spatial resolution of 
the 53 Mpc cases also yields somewhat higher mean density of the 
ionized regions compared to the corresponding 163 Mpc box cases. 



4.3 The effects of the source model: photon production 
efficiencies and minimum source mass 

In Fig.[9](right panels) we show the corresponding reionization his- 
tory results when the source models are varied. We also replotted 
one of our fiducial cases, S 1 , for facilitating direct comparison with 
the self-regulated cases. All reionization histories remain largely 
monotonic throughout the evolution, which therefore is a fairly ro- 
bust feature, independent of the particular ionizing source proper- 
ties assigned. However, a wide range of overlap epochs - from as 
early as z — 12.9 in the no suppression case S3 to z — 8.3 in large- 
source-only case S5, and a wide range of slopes of the reionization 
history evolution are observed. 

Our fiducial case, SI, has the most extended reionization his- 
tory of all, which starts with the formation of the first 10 8 Mq halos 
at z ~ 26 and reaches overlap at z = 8.9. In comparison when all 
sources have the same efficiencies, but none are ever suppressed 
(case S3) the reionization history is very steep, roughly exponen- 
tial, tracking the exponential rise of the collapsed fraction in ha- 
los (cf. Figure[3]l. The no-suppression, low-efficiency case S4 also 
produces a very extended history since it also starts with the for- 
mation of the first 10 8 Mq halos and, by design, reaches overlap at 
roughly the same time as SI. However, the sources are necessar- 
ily much weaker in that case compared to SI and S3, and there- 
fore the ionized fraction starts much lower compared to the fiducial 
simulation and only catches up with the self-regulated case at late 
times (x m J> 0.25). Finally, in case S5 only the massive sources 
are active, and therefore the reionization starts late, at z ~ 18, 
but x m rises exponentially, in proportion of the collapsed fraction 
in those massive halos, reaching (again by design) overlap at the 
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Figure 10. (bottom panels) Number of ionizing photons emitted by all sources (i.e. if there were no suppression; thin lines) and all active sources (thick lines) 
in the computational volume per timestep and (top panels) cumulative number of photons per total gas atom released into the IGM. Notation is the same as in 
Fig. [9] Vertical lines with same colors and linetypes mark the overlap redshift in each case. 



same time, z = 8.3 as the fiducial case. The lack of low-mass ha- 
los therefore delays reionization considerably and naturally yields 
much lower intergated electron scattering optical depth (0.07 1 com- 
pared to 0.084 for SI, and 0.078 for S4). 

The mean overdensity of the ionized regions, x m /x v (Fig. [9] 
top right panel) is above unity for all cases and at all times, demon- 
strating the robustness of the inside- out nature of reio nization, in 
agreement with our original findings diliev et alJuOOoi) . Compared 
to our fiducial simulation, SI, the ratio x m /x v is significantly 
higher for cases S4 and S5. This is due to the fact that the H II 
regions are more tightly correlated with the density peaks in those 
cases, because the number of sources, which form at the density 
peaks, rises exponentially in these cases. Finally, simulation S3 
show an intermediate behaviour, similar to the fiducial case, SI, 
but with somewhat faster decrease of the mean overdensity of the 
ionized regions, due to higher ionization of the low-density regions 
in this case. 

The corresponding evolution curves of the cumulative num- 
ber of ionizing photons emitted within each simulation are shown 
in Fig. [TO] Indicated are also the number of photons which would 
have been emitted if no self-regulation has taken place (thin lines) 
and the overlap redshifts (vertical marks). Starting with our self- 
regulated cases (Fig. [TO] left) we see that the low-mass source sup- 
pression does not have a significant impact until 2 ~ 20 — 22, but 
after that has a great impact, reducing the overall number of emit- 
ted photons by up to a factor of ~ 30 for highly-efficient low-mass 
sources and about a factor of 10 for less efficient ones. The number 
of ionizing photons emitted per timestep also becomes variable and 
non-monotonic function of redshift due to the complex interplay of 
suppression and new source formation during the self-regulation 
process. Eventually, by overlap (indicated for each case by the ver- 
tical lines) all low-mass sources are suppressed and the number of 
photons continues to rise smoothly as ever more high-mass sources 
form. The simulation boxsize makes little difference in the ioniz- 
ing photon production, apart from modest variations due to the dif- 



ferent random realization in each case. In terms of the cumulative 
numbers of emitted photons per atom (top panel), by redshift 2 = 6 
up to 10 photons per atom are produced in the efficient-sources case 
and up to 2 photons per atom in the low-efficiency cases. However, 
at their respective overlap redshifts approximately the same number 
are produced, about 1.5, i.e. on average only about one recombina- 
tion per every 2 atoms occurs during the evolution. Therefore, re- 
combinations are relatively unimportant in these runs. The reason 
for this is that much of the density fluctuations are at very small 
scales, well below our radiative transfer grid resolution. This ad- 
ditional small-scale power can be added as sub-grid clumping of 
the gas, calculated based on much higher resolution simulations. 
We will consider the effects of sub-grid clumping in a companion 
paper (Koda et al., in prep.). 

Turning our attention to the set of cases with different source 
efficiency models (Figure QJJ right), we first note that the three 
models which by construction have very similar overlap epochs 
(our fiducial case, SI, and cases S4 and S5) also have almost iden- 
tical photon production numbers at overlap. All three reach this 
point in a very different manner, however. In run S4 all sources 
are quite weak, but no sources are ever suppressed, and the emis- 
sivity per timestep reaches the fiducial case once most low-mass 
sources are suppressed in the latter case. In case S5 only the mas- 
sive sources are present and consequently its emissivity lags signif- 
icantly at early times until eventually the exponential rise of those 
sources allows it to join the other two cases at 2 <. 11. We also note 
that after overlap the photon emissivity in the no-suppression case 
S4 lags behind the others because its high-mass sources are very 
inefficient and the collapsed fraction in low-mass sources by this 
point does not rise as fast as the one for the high-mass sources. Fi- 
nally, the number of photons produced in the no suppression, high 
source efficiency case, S3, simply follows the total collapsed frac- 
tion in all sources and therefore rises almost exponentially, roughly 
parallel to the curve for S4, eventually surpassing 10 photons per 
atom before z — 10. However, at its own (very early) overlap epoch 
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Figure 11. Spatial slices of the ionized and neutral gas density from our radiative transfer simulations with boxsize 163 Mpc: (a)(left) LI (b)(middle) L2, 
and (c)(right) L3, all at box-averaged ionized fraction by mass of x m ~ 0.50. Shown are the density field (green) overlayed with the ionized fraction 
(red/orange/yellow) and the cells containing active sources (dark/blue). 




Figure 12. Spatial slices of the ionized and neutral gas density from our radiative transfer simulations with boxsize 53 Mpc, all at box-averaged ionized fraction 
by mass x m ~ 0.50. Shown are the density field (green) overlayed with the ionized fraction (red/orange/yellow) and the cells containing sources (dark/blue). 
Shown are (left to right and top to bottom) cases SI, S2, S3, S4, and S5. 



even this case produces the same number of photons as the others, 
about two per IGM atom. 

In Figures [7TI and [72l we show slices through our simulation 
volume showing the geometry of the H II regions at x m ~ 0.5, 
overlayed on the corresponding density field for all our simula- 
tions. We also mark the cells containing active sources (blue/dark). 
Comparing first the large, 163 box cases (Fig. II lb , we note that, as 
could be expected in the high-efficiency fiducial case LI there are 



many more active sources in the low-efficiency one, L2. In contrast, 
in simulation L3 there are many fewer sources dure to its higher 
mass cutoff, which only leaves the high-mass, rare sources present. 
The large-scale structures are quite similar in size and shape in 
the two self-regulated simulations, but the fiducial case LI yields 
much more small-scale ionized patches even though it reaches 
half-ionized state noticeably earlier, at which point there are many 
fewer, and more clustered, sources. The reason for this apparently 
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counter-intuitive behaviour is that the much weaker sources in case 
L2 have difficulties fully ionizing their own cells (which, at 445 
kpc/h linear size, are relatively large), and therefore large number 
of sources are needed to produce a sizeable fully-ionized patch. In 
contrast, the much more efficient sources in our fiducial case LI 
easily ionize their own cell, resulting in many small-size H II re- 
gions, instead of the more scattered, partially-ionized cells in L2. 
On the other hand, in the non-self-regulated case L3 we find many 
fewer, la rger ionized regio ns, in agreement with our previous re- 
sults in dlliev et al.ll2008af) . The large-scale structures have some 
similarities to the ones found in the self-regulated cases, as could 
be expected given that all simulations share the underlying large- 
scale cosmic structures. However, the ionized regions are in a more 
advanced stage where they start merging together, and there is far 
less small-scale structure due to the absence of low-mass, weaker 
sources. In that case there are also no partially-ionized regions, 
since there is no low-mass source suppression Jeans-mass filtering 
(i.e. sources do not die), and all sources are sufficiently luminous 
to completely ionize their own local volume. 

The corresponding images from our small-box simulations at 
the same ionized fraction of x m ~ 0.5, are shown in Fig.[T2] Once 
again, the large-scale structures, which tend to strongly corellate 
with the underlying distribution of density and clustered halos, are 
generally quite similar. In contrast, ionized patches produced by 
multiple, less biased sources whose distribution does follow the 
knots and filaments of the Cosmic Web are much more irregularly- 
shaped. There are significant differences in the smaller-scale struc- 
tures among the range of simulations. The self-regulated cases, 
SI, S2 and to a lesser extend case S4 have the most small-scale 
structure, including both small H II regions and rough, irregularly- 
shaped large H II region boundaries. In contrast, S5, which does 
not include the low-mass sources and S3, which includes efficient, 
unsuppressible low-mass sources both yield many fewer ionized 
patches with smoother boundaries, which reflects the rarity and 
highly clustered nature of their active sources. We have presented 
a more detailed discussion of the H II region geometry, size distri- 
bution and topological characteristics in a recent co mpanion paper 
based on a subset of the current suite of simulations dFriedrich et al.l 

I201I . 



5 OBSERVATIONAL SIGNATURES 

We now turn our attention to the reionization observables and 
specifically how are they related to the assumed source populations 
and their efficiencies. A better understaning of these dependencies 
should allow us in turn to use the observational data to constrain 
the properties of the reionization sources. Our main focus will be 
on the redshifted 21 -cm signatures, although we also briefly discuss 
the photoionization rates in the IGM, related to the measurements 
of the Gunn-Peterson effect and the gas temperature. 



5.1 Photoionization rates 

An important, if indirect, observable signature of reionization is the 
mean photoionization rate in the IGM. At present this quantity has 
only been measured for the post-reionization IGM at z < 6, derived 
based on the small residual neutral fraction and its corresponding 
Ly-a optical depth. It therefore typically characterizes only the fi- 
nal EoR stages, around and after overlap. 

The redshift evolution of the mean photoionization rates, F, 
averaged over our simulation volume for our large-box simulations 
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Figure 13. Evolution of the mean mass-weighted (thin lines) and volume- 
weighted (thick lines) photoionization rates in our computational volume 
for simulations LI (red), L2 (blue) and L3 (black). 
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Figure 14. Photoionization rate PDF's for our fiducial case LI for epochs 
when the ionized fraction by mass is x m = 0.1 (red), x m = 0.5 (blue), 
^m = 0.9 (green) plotted in linear (top) and log (bottom) scales. 



L1-L3 are shown in Figure Q~3] In overall curve shape and timing 
the evolution roughly mirrors the reionization histories for these 
three cases. This is natural, since the mean photoionization rate is 
average of the fraction from the ionized regions, where the F values 
are high and fairly uniform at ~ 10 -12 s _1 (see also Fig.ll4b, and 
the neutral regions, where r ~ 0. The mean photionization rate in 
our fiducial case LI initially rise roughly exponentially, until self- 
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regulation becomes wide-spread at ~ 16, at which point the stall 
at around 10~ 14 until the high-mass sources become sufficiently 
abundant to dominate the evolution, which occurs around z ~ 11. 
At this time the mean rates resume their steep rise, reaching a peak 
of r ~few x 10 -12 by overlap. The situation is somewhat differ- 
ent for the low photon production efficiency case, L2. The mean 
r values are much lower in this case, by 2-4 order of magnitude, 
and are mostly rising monotonically throughout the evolution, as 
the Jeans-mass suppression effects are milder in this scenario. The 
peak value reached in this case is about 1CP 12 , in rough agreement 
with the measured one at z ~ 6. However, we should note that any 
direct comparisons to observationally-derived values are at best ap- 
proximate since our simulations at present do not take into account 
the Lyman-limit systems (LLS), which are likely to limit the growth 
of the mean free path of the photons and thus limit T, as well. Be- 
fore overlap the mean free path is dictated by the remaining neutral 
regions and the (still fairly high) residual neutral gas fraction within 
the ionized regions, and therefore the LLS are unimportant and do 
not affect our simulation results. The same is probably not true after 
overlap and we will study the effects of LLS in future work. For our 
current purposes the lack of LLS means that we cannot yet make 
a firm conclusion that the low efficiency, late-overlap case, L2, fits 
the observations better than case LI. 

The mean photoionization rate for the non-self-regulated case 
L3 is intermediate between LI and L2. It starts from very low 
values, ~ 10" 19 , when there are still only a very few high-mass 
sources, around z ~ 17, but then rise sharply, roughly expo- 
nentially, and converges (by construction, since efficiencies were 
picked so they overlap at the same time) to the values for LI at 
later times. 

The mass-weighted photoionization rates (thin lines) are sig- 
nificantly higher, by factors of up to 2-3 than the volume-weighted 
ones (thick lines) at all times and for all simulations. This is easy 
to understand given the inside-out nature of reionization, whereby 
the ionizing sources are found in dense regions, which pushes the 
mass-weighted means higher. Such large differences are interest- 
ing, however, since they can possibly skew the observationally- 
derived values. Probes of the mean, low-density IGM will there- 
fore yield considerably lower values for Y than any measurements 
which are more sensitive to denser regions, e.g. around sources. 

Several illustrative PDF's (at cell size, here 445 kpc/h) of the 



photoionization rates for our fiducial simulation LI are shown in 
Figure [74l Plotted are the PDF's at early (x m — 0.1), intermediate 
(Xm — 0.5) and late (x m — 0.9) stages of the evolution. These can 
be compar ed to the no low- mass sources data we presented pre- 
viously in llliev et alJ l2008b). At all times the PDFs show a char- 
acteristic, three-peaked profile. The rightmost peak, at T_i2 ~ 1, 
is formed by the cells inside the H II regions, while the other two 
peaks correspond to partially-ionized cells, predominantly at the 
expanding I-fron ts and relic (i.e. recombining) H II regions. As we 
have shown in (Iliev et al. 2008b), the ionization state is close to or 
at equilibrium deep inside the ionized regions, but far from equilib- 
rium at the I-fronts. This holds true for the current simulations, as 
well. However, compared to our previous simulations the low-mass 
source suppression in the current runs yields a significant fraction 
of volume in relic H II regions and partially-ionized cells and thus 
higher peaks at lower T values than was observed in the simulations 
without suppression. 

In Figure Q3] we show scatter plots and the corresponding 
contour levels of the local photoionization rates, T vs. density in 
units of the mean, 1 + 5 = p C cii/p for our fiducial case, LI and 
Xm = 0.1, 0.5 and 0.9. Overall, there is a clear positive correlation 
between the density and the photoionization rate. This could be ex- 
pected, given that sources, around which the photorates peak form 
preferentially in high-density regions. However, the relationship 
between the two is complex, the correlation is weak and the scat- 
ter significant. Similarly to the PDF's discussed above, three peaks 
are observed. The high-r peak (F ~ 10~ 12 s _1 ) consists of the 
ionized cells, which are typically denser than average. The middle 
peak, at F ~ 10 -15 s _1 , corresponds to I-fronts and other partially- 
ionized regions, while the cells with still lower values (r ~ 0) 
correspond to the still-neutral regions. At early times (x m — 0.1) 
the majority of cells is either neutral or partially-ionized and the 
correlation with the local density is very weak. When the process 
advances (x m = 0.5) a large population of fully-ionized, high- 
r cells develops and within the H II regions the photoionization 
rate is fairly well correlated with the density, albeit still with a 
large scatter. On the other hand, for the photoionization rates in 
the partially-ionized regions shows essentially no correlation with 
the density. At late times (x m = 0.9) these trends become even 
more pronounced and a quite tight correlation develops for the 
highest-density regions. The overall behaviour is consistent with 
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what w e previously obse rved when no source suppression were 
present jlliev et alj2008a) . albeit with some minor quantitative dif- 
ferences. 



5.2 Redshifted 21-cm 

The differential brightness temperature of the redshifted 21-cm 
emission with respect to the CMB is determined by the density of 
neutral hydrogen, pui, and its spin temperature, Ts and is given by 

Ts — Tcmb 



Ts - Tcmb 3X 3 A 10 T,n H i{z) 



1 + z 
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I IFieldl ll959). where z is the redshift, Tcmb is the temperature of 
the CMB radiation at that redshift, r is the corresponding 21-cm 
optical depth, assumed to be small when writing equation [3] Ao = 
21.16 cm is the rest-frame wavelength of the line, Aio — 2.85 x 
10 -15 s _1 is the Einstein A-coefficient, T* = 0.068 K corresponds 
to the energy difference between the two levels, 1 + 8 — nui /{uh} 
is the mean number density of neutral hydrogen in units of the mean 
number density of hydrogen at redshift z, 

f2bPcrit,0 



(n H )(z) 



iiH"m v 



-(i + *)' 



= 1.909xl0- 7 cm- 3 (^)(l + z) 



(4) 



with /in = 1.32 the corresponding mean molecular weight (as- 
suming 24% He abundance), and H(z) is the redshift-dependent 
Hubble constant, 



h(z) = # [fim(i + zf + n k (i + zf + n A ] 1/2 

= H E(z)^H nll 2 (l + zf /2 , 



(5) 



where Ho is its value at present, and the last approximation is valid 
for z 2> 1. Throughout this work we assume that Ts 3> Tcmb i.e. 
that all of the neutral IGM gas is Ly-a-pumped by the background 
of UV below 13.6 eV from early sources and heated well above the 
CMB temperature (due to e.g. a small amount of X-ray heating), 
and thus the 21-cm line is seen in emission. These assumptions are 
generally well-justified, except possibly at the earliest times (see 
e.g. lFurlanetto et alj2006bt and references therein). 



5.2.1 Evolution of the patchiness 

In Figure [76l we show space-redshift (space-frequency) slices cut 
through the simulated image cube as a radio array would see it 
(ignoring foregrounds). The spatial dimension is on the vertical, 
where we have duplicated the computational volume for visualiza- 
tion purposes, while the redshift/frequency is along the horizontal. 
Images are of the 21-cm emission differential brightness tempera- 
ture signal extracted from our fiducial simulation LI, continuously- 
interpolated in redshift/frequency including redshift-space distor- 
tions due to peculiar velocities. The volume is cut at an oblique 
angle in order to minimize artificial repetition of structures along 
any line of sight. The top and middle panel show images (in log, 
which shows better the residual H I fraction in the ionized regions 
and linear scale, which shows better the neutral structures) at the 
full simulation resolution, which is much higher than what current 



experiments will achieve given the sensitivity constraints. The bot- 
tom panel shows same data, but smoothed with a Gaussian beam 
and an integrated bandwidth which both roughly correspond to the 
values adopted in the LOFAR EoR experiment. To mimic the fact 
that an interferometer such as LOFAR is insensitive to the global 
signal, the mean signal at every frequency slice has been subtracted. 
At high redshift, here z > 11 all H II regions are small and 
largely isolated. Smoothing the data to the LOFAR resolution (the 
ones for MWA and GMRT are even lower) renders such small 
structures undetectable. One needs at least ~ 1' or better resolu- 
tion for potentially observing them, making this regime a potential 
target for future, more sensitive experiments e.g. SKA. However, 
at intermediate redshifts (here z ~ 10) the ionized regions quickly 
grow by merging and remain clearly visible also after beam- and 
bandwidth smoothing. Even though some detail is lost, the large- 
scale structure of the ionization field remains visible all the way to 
the overlap epoch, here z = 8.4. As was noted above, compared 
to the simulations with no self-regulation (e.g. Iliev et al. 2008b), 
the suppression of low-mass sources introduces much more small- 
scale structure and many, mostly small partially-ionized and relic 
H II regions. However, the smoothing to the radio array resolu- 
tion largely eliminates this fine-scale structure and the result is, at 
least visually, not dramatically different from the case with no self- 
regulation. The minimum and maximum values of the differential 
temperature are also similar. We consider more quantitative mea- 
sures of the 21-cm signal next. 

5.2.2 21-cm background: mean and rms 

The evolution of the mean differential brightness temperature and 
its rms fluctuations for our fiducial case LI and L3, which cor- 
responds to our p revious simulations with no self-regulation in 
diliev et alj |2008b) are shown in Figure [17] The presence of low- 
mass sources and Jeans mass filtering yields initially a steeper de- 
cline of the mean 21-cm emission starting from v ~ 80 MHz (z ~ 
17), at which point the low-mass sources start forming in larger 
numbers (becoming v ~ 3 halos, cf. Figure[3}, while the high-mass 
sources are still very rare. At v ~ 130 MHz (z ~ 10) the high- 
mass sources in turn become 3 — a halos, i.e. relatively more com- 
mon and the mean STt evolution for case L3 steepens, eventially 
reaching the same overlap epoch (by construction). In terms of de- 
tectability in experiments looking for looking f or rapid changes 
in the 21cm signal as the Universe reionizes JShaver et all 1 19991 : 
iBowman & RogerjuOlOT) . this behaviour means that the case of 
self-regulation is even more difficult to detect than the one without. 
While without self-regulation the global signal drops fast by about 
25 mK between 130 and 150 MHz, with it the drop at the higher 
frequencies is more gradual. The decrease with self-regulation is 
somewhat steeper at lower frequencies, v — 80 — 120 MHz, but it 
is still fairly gradual and more difficult to detect. 

Comparing the rms fluctuations averaged over LOFAR-like 
beam and bandwidth (Figure[l7] right) we see that the overall evo- 
lution follows similar paths in both cases. Early-on very little of the 
gas is ionized and the fluctuations therefore simply follow the den- 
sity ones. Only when a significant ionized fraction develops do the 
fluctuations depart from the underlying density. For simulation L3 
this occurs fairly late, at v > 110 MHz, compared to much earlier, 
v > 80 MHz, for the fiducial simulation. At this point the rms fluc- 
tuations slightly dip, as the highest density peaks are ionized, which 
diminishes the mean STt but does not boost the fluctuations since 
the H II regions are still smaller than the smoothing size. As the H II 
regions grow, the fluctuations increase again, reaching a peak be- 
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Figure 18. The evolution of the mean 21-cm background for our fiducial cases (left) and varying the source model (right). All cases are labelled by color and 
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fore the signal dips again as the IGM becomes highly ionized. The 
peak position (at 142 MHz) remains the same in both cases and is 
thus not affected by self-regulation. The rms fluctuations are lower 
with self-regulation, by about 1/3 at the peak. The reason for this is 
the lower mean differential brightness temperature in that case, as 
can be seen in the bottom panel. When the fluctuations are normal- 
ized by the mean they become indentical in the two cases once they 
surpass the density fluctuations (y > 127 MHz). Before that point 
the fluctuations in case L3 closely follow the density ones, while 



in case LI they are lower because the highest-density peaks have 
already been ionized. 

In Figure [T8] we show the evolution of the mean redshifted 
21-cm differential brightness temperature for high vs. low source 
efficiencies and different box sizes (left) and for varying source 
models (right). The lower photon efficiencies (simulation L2) pre- 
dictably yield a more gradual transition of the global IGM from 
neutral to ionized state compared to our fiducial case LI. E.g. the 
evolution from 25 mK to ~ mK occurs over ~ 50 MHz, from 
130 to 180 MHz. Such an evolution would make detection of the 
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'global step' even more difficult compared to the fiducial case. The 
computational volume adopted for the simulation makes very lit- 
tle difference to the predicted mean 21 -cm signal, demonstrating 
again that 37 Mpc/h box is sufficiently large to faithfully represent 
the mean reionization history. Most cases with varying UV-source 
models (Figure Q~8] right) yield mean 21 -cm histories which are 
quite similar to each other, a consequence of their analogous reion- 
ization histories. The only noticeable differences are at intermedi- 
ate frequencies, between 90 and 130 MHz, where case S4 gives 
higher STt by up to 5 mK. The only significantly different evo- 
lutions are provided by cases S5 and S3. Those scenarios exhibit 
sharper 21 -cm step due to their faster, exponential rise of the ion- 
izing photon emissivity due to the weak or no suppression in those 
cases. 

The evolution of the rms 21-cm emission fluctuations for 
LOFAR-like beam and bandwidths corresponding to the same sets 
of simulations as in Figure [T8] are shown in Figures [T9] (high vs. 
low ionizing efficiencies and varying boxsize) and[20] (different UV 
source models). We show the same data vs. observed frequency 
(left) and ionized fraction (right). The latter takes away the reion- 
ization timing and allows comparison at the same stages of each 
reionization history regardless of when they actually occur in time. 

In all cases the rms evolution roughly follows the same path, 
with an initial rise tracking the underlying density fluctuations 
when the IGM is mostly neutral, with a subsequent decrease when 
the first H II regions appear followed by a second, higher peak of 
the fluctuations at later times when the initially small ionizing re- 
gions grow, overlap locally and as a result match better the inter- 
ferometer beam and bandwidth resolution, and a final decline when 
most hydrogen is ionized. However, despite this recurring pattern, 
there are significant, interesting, and often instructive differences 
among the models. Varying the ionizing photon emission efficien- 
cies primarily changes the timing of the peak of the fluctuations 
(Figure[T9l left), from 142 MHz (z — 9) for the fiducial case LI to 
172 MHz (z = 7.24) for the low-efficiency case L2. However, as 



seen in Figure [19] (right) this shift to later times is fully explained 
by the delayed reionization in the low-efficiency model and both 
curves peak at mass- weighted ionized fraction x m ~ 0.7. We note 
that the latter value is dependent, apart from the reionization pa- 
rameters, also on the beam and bandwidth considered, as the peak 
is reached when the typical H II region size is best matched to the 
radio array resolution, therefore the peak occurs earlier in the reion- 
ization history for higher resolution and later for lower one. The 
fluctuation dip due to the earliest H II regions also occurs at the 
same point of the reionization history (x m ~ 0.2 — 0.25) in both 
cases, but the lowest rms values reached differ significantly at 1.5 
mK for the fiducial model vs. 2.3 mK for the low-efficiency one. 
The reason for this is that in the former case the bottom occurs ear- 
lier, when the sources responsible are rarer, more biased and there- 
fore ionize the highest density peaks, which results in a larger de- 
crease in the fluctuations. The simulation volume (and, correspond- 
ingly, radiative transfer grid resolution) has moderate, but apprecia- 
ble effect on the 21-cm fluctuations in our fiducial case. The peak 
height is decreased by 7%, from 6.1 mK to 5.7 mK, but it is also 
shifted to earlier time/lower frequency (to 138 MHz). Interestingly, 
while the dip of the fluctuations is also shifted to lower frequency 
for the smaller box simulation, the lowest rms value is in fact 
higher. At first sight this appears counter-intuitive, since naively 
we might expect that the higher grid resolution in the smaller vol- 
ume to yield a larger rms decrease (since the density peaks where 
the first sources form are resolved better in this case). What actu- 
ally occurs here is more complicated, however. Statistically, there 
are fewer (and lower) high-density peaks in the smaller volume, 
which diminishes the effect of the very first sources on the fluctu- 
ations. Furthermore, the very first sources form later in the smaller 
box, again due to its much smaller volume, which means that the 
21-cm rms fluctuations track the density ones for somewhat longer. 
On the other hand the effects of box size and resolution for the low- 
efficiency cases are small and manifest themselves solely through 
the higher underlying density fluctuations in the small-box, high- 
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Figure 20. Same as Fig. 1 191 but for varying UV source models: SI (green, short-dashed), S4 (cyan, long dashed), S5 (light red, dot-short dashed), and S3 
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Figure 21. The effect of self-regulation on the 21-cm differential brightness temperature fluctuation power spectra. Shown are the epochs at which the ionized 
fractions are (left) x m = 0.2, (middle) x m = 0.5 and (right) x m = 0.75 for our fiducial self-regulated case, LI (blue, solid) and the corresponding 
non-selfregulated case with same overlap epoch, L3 (black, long-dashed). 



resolution simulation. Both the peak and dip reach the same value 
for the two cases and occur at the same frequencies. 

On the other hand, variations of the ionizing source model, 
yield a wider variety of 21-cm rms evolutions (Figure l20t. Inter- 
estingly, all models exhibit the basic evolution features seen in our 
fiducial simulations - the initial dip of the rms value when the first 
H II regions appear, followed by a (relatively narrower) peak at 
later times when the process is sufficiently advanced for the typ- 
ical patch size to roughly match the radio beam and bandwidth 
(though we note that more extreme, and unrealistic, source mod- 
els can produce rms fluctuations with a very different shape, see 
Appendix lAt. There are significant variations in the details of the 
evolution, however. The simulations with no low-mass source sup- 
pression but same overlap as in our fiducial case (S4 and S5) yield 
rms fluctuations peaks which are at roughly the same frequency 



as the fiducial case (v ~ 140 MHz, more specifically), but the 
(ST£) 1/2 peak values are up to 50% higher (~ 7 - 8 mK). The 
early reionization, no suppression case, S3, gives a still higher peak 
value, reaching almost 10 mK, and a narrower peak. 

A different way to consider the same data is to plot the differ- 
ential brightness temperature evolution in terms of its own reion- 
ization history, i.e. against the mass-weighted ionized fraction, x m , 
(Figure [20] left), which removes the dependence on the absolute 
timing of reionization. There is only a modest variation in the reion- 
ization stage (i.e. ionized fraction, x m ) at which the rms peak is 
reached, which ranges x m ~ 0.6 — 0.7 (i.e. relatively late in the 
reionization history). However, as we noted above, the peak value 
itself varies significantly between the simulations. It is highest for 
the high-efficiency, no suppression model, S2 ((ST 2 ) 1 ' 2 = 9.4 mK 
compared to 3.8 mK for our fiducial case SI). Lower efficiency and 
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no suppression (S4) results in a significantly lower peak (6.9 mK), 
while when only massive halos are present (L5) the peak is again 
quite high {(5Tb) 1 = 8.0 mK). Therefore, in general more 
abrupt reionization scenarios (S3, S5) result in higher fluctuations 
at the peak, while more extended ones (due to self-regulation or 
lower efficiencies) give lower rms peak values. 



5.2.3 21-cm background fluctuations: power spectra 

We now turn our attention to the (3D) power spectra of the 21-cm 
emission derived from our simulations. We construct the brightness 
temperature datacube in the redshift space using what we term the 
PPM-RRM scheme, as follows. We first develop an adaptive-kernel, 
SPH-like approach to compute the bulk-flow velocity of the IGM 
at any position, directly from N-body particle data. We paint the 
particle mass by the hydrogenic neutral fraction of the RT grid that 
the particle resides. Then N-body particles are Doppler-shifted to 
their apparent locations by LOS bulk-flow velocity, new smooth- 
ing kernel lengths are computed using the new particle positions in 
redshift-space, and halo-excluded particle data (i.e., H I mass) are 
again smoothed onto a regular, redshift-space grid at RT grid reso- 
lution. Then we compute the redshift-space HI density fluctuation, 
and 21-cm brightness temperature measured in redshift space by 



sn(s) = m( Zcos ) [i + s; H1 (s) 



(6) 



where 5Tb was defined in equation [3] We compute the redshift- 
space power spectrum using FFT. We refer the reader to lMao et al.l 
1 1201 II) for a detailed discussion of this methodology. 

In Figure[2T]we compare the 21-cm power spectra for models 
L3 (equivalent to our previous simulations with no self-regulation) 
and our fiducial self-regulated case LI for several representative 
stages of reionization. In both cases the 21-cm power spectra are 
initially close to a power law, with no characteristic scale, and with 
only the L3 power flattening out at small scales due to its lack of 
very small-scale structures and smooth, large H II regions. Overall, 
there is always less power in our fiducial case, by factor of 50% (in 
Afc), or more. As ionized regions continue to expand a character- 
istic scale starts to emerge, which for these particular simulations 
is around wavenumbers k ~ 0.2 — 0.8 h/Mpc. Interestingly, this 
feature shows up at the same scales regardless of the presence of 
low-mass sources, indicating that the characteristic H II region size 
is caused by the clustering of the high-mass, unsupressible sources. 

On the other hand, the assumed ionizing source efficiencies 
have at most only a minor effect on the power spectra once the self- 
regulation is included (Figure [22l. The characteristic H II region 
scale is the same in the two cases and arises at the same point in the 
reionization history. The modest differences in the power spectra at 
the early stages of reionization arise due to the preferential ioniza- 
tion of the high density peaks (where the first sources form). At the 
same ionized fraction x m there are many more sources in the low- 
efficiency case, forming very small H II regions, vs. fewer, larger 
ones in the high-efficiency case. As a result, the low-efficiency 
model yields less power at very small scales (k > 4), but a little 
more power at intermediate scales (k — 0.2 — 2). At the middle 
and late stages of reionization the two power spectra at the same 
x m are largely identical, with only small differences due to the 
different amount of small-scale structure. The simulation volume 
utilized also has little effect on the derived power spectra at early 
times, essentially just shifting the range of k over which the results 
obtained are reliable. However, as the H II regions grow at interme- 
diate and late stages of reionization, their sizes become comparable 



to the simulation volume and as a result the 37 /h Mpc volume can- 
not represent the bubble size distribution correctly and the calcu- 
lated power spectra completely miss the characteristic H II region 
scale. On the other hand, the smaller volumes yield more fluctua- 
tion power at small scales, due to their superior spatial grid reso- 
lution. These results argue for a strong caution when using small 
(sub-100//i Mpc) boxes for predicting any EoR observables at late 
times. 

Turning our attention to the varying source models (Fig- 
ure 123b . several trends become clear. The lack of Jeans mass fil- 
tering (case S3 vs. the fiducial SI) results in much more power at 
large scales (k < 5 h/Mpc early, k < 2 h/Mpc at late times), but 
considerably less power on small scales, consequence of the large 
H II regions with smooth boundaries in case S3 produced by its 
luminous and strongly clustered sources. This also results in very 
flat power spectra for S3, with roughly constant power at all scales 
at intermediate and late times (x m = 0.5 and 0.9). On the other 
hand, if the low-mass sources are not present at all (case S5 vs. S3) 
there is more power on all scales during the early stages of reion- 
ization (x m — 0.1). However, at intermediate and late stages of 
reionization the power spectra for the large-source only case S5 re- 
main steeper, with considerably more power at small scales, which 
results in cases S3 and S5 power spectra crossing at k ~ 2 h/Mpc. 
The reason for this somewhat counterintuitive behaviour is that the 
same reionization stage is reached in S5 much later than in S3, 
by which time there are many more and less clustered high-mass 
sources, which in turn results in more power at small scales. Fi- 
nally, the low-efficiency case with no self-regulation (S4) yields 
very similar power spectra to our fiducial simulation S 1 throughout 
the evolution, with only slightly more power at large scales dur- 
ing the early and intermediate stages, and slightly less small-scale 
power at the late stages. This suggests that the SI and S4 scenar- 
ios might be difficult to distinguish solely through power spectra 
measurements. 



5.2.4 21-cm background fluctuations: PDFs and 
non-Gaussianity 

The 21-cm power spectra would fully characterize the emission 
field if the differential brightness distribution were purely Gaus- 
sian. However, generally that is n ot the case during reionization, as 
we have shown in previou s work JMellema et al. 2006b; Ilie v et al.l 
l2008al ; iHarker et al.ll2009h . The probability distribution functions 
(PDFs) of 5Tb could be significantly non-Gaussian, p articularly 
at the later stages of reionization JMellema et al. 2006b) and their 
measured skewness ca n be used to discrim inate between different 
reionization scenarios Marker et alj 120091) . The PDF's and their 
evolution could also be used to derive the reionization history of 
the IGM dlchikawa et al.ll2010l : lGluscevic & Barkanj|2010l) . 

The 21-cm cell-by-cell PDFs with and without the presence 
of low-mass, suppressible sources at three representative stages of 
reionization (x m — 0.1, 0.5 and 0.9) are shown in Figure[26] Early 
on (x m = 0.1) the distributions are mostly following the underly- 
ing density field, and as a consequence are mostly Gaussian. There 
is a non-Gaussian tail for high differential brightness temperatures 
due to density nonlinearities. Reionization itself introduces some 
non-Gaussianity at low 5Tb as the first H II regions form around 
the highest density peaks, moving the corresponding cells into the 
extreme left of the distributions (i.e. holes in the neutral gas dis- 
tribution, with STt ~ 0). This slightly skews the distribution to- 
wards below-average (i.e. negative in 5Tb — 5Tb) temperature val- 
ues since the low-density regions remain more neutral, on average. 
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Figure 22. 21-cm differential brightness temperature fluctuation power spectra. Shown are the epochs at which the ionized fractions are (left) x m = 0.2, 
(middle) x m = 0.5 and (right) x m = 0.75. All cases are labelled by color and line-type, as follows: LI (blue, solid), SI (green, short-dashed), L2 (magenta, 
long-dashed), and S2 (red, dotted). 
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Figure 23. 21-cm differential brightness temperature fluctuation power spectra for varying source models. Shown are the epochs at which the ionized fractions 
are (left) x m = 0.1, (middle) x m = 0.5 and (right) x m = 0.9. All cases are labelled by color and line-type, as follows: SI (green, short-dashed), S4 (cyan, 
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Figure 25. Same as in Fig. 1241 but for boxcar smoothing of 5 h 1 Mpc. 
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Figure 26. The effect of the source efficiencies (high- vs. low efficiency) and box size on the PDF distribution of the 21-cm signal. Shown are the epochs at 
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Figure 27. Same as in Fig. 1261 but for boxcar smoothing of 5 h 1 Mpc. 
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Self-regulation mitigates those trends somewhat, as the low-mass 
sources are less clustered and more uniformly distributed through- 
out the volume, rather than being only at the highest density peaks 
(which are strongly clustered, as a consequence of the Gaussian 
density field statistics, see Figure[6}. For the same reasons the PDF 
with self-regulation is also slightly less wide than without, as evi- 
denced by the Gaussian distributions with the same mean and width 
as the actual PDFs. During the later stages of reionization the distri- 
bution becomes ever more non-Gaussian, with the most prominent 
feature due to the ionized regions (5Tb — 5Tb < 0). The remaining 
neutral regions are a mixture of voids (low, but positive 5Tb — 5Tb) 
and a few remaining higher-density regions (e.g. filaments) which 
form the 21-cm bright non-Gaussian tail. Once again, both features 
are much reduced when low-mass source suppression is taken into 
account due to the weaker clustering of such sources. The width of 
the PDFs decreases over time, and does so faster when the low mass 
sources are not present. At late times the two distributions have al- 
most the same means and widths, although the actual distributions 
remain significantly different. 

When the same PDFs are smoothed with a 5 Mpc/i -1 window 
(roughly similar in size to e.g. the LOFAR beam, albeit here we 
use a different window shape for simplicity of the calculations) the 
results become notably different (Figure [2Tb. The smoothed PDFs 
at early times (x m = 0.1) become significantly more Gaussian, al- 
though some residual non-Gaussian tails remain at both high and 
low 5Tb. As could be expected, the smoothed distributions also be- 
come much less wide compared to the unsmoothed ones, since the 
smoothing window averages the values, flattening the highest peaks 
and deepest valleys of the distributions. 

Interestingly, at the middle and late stages of reionization 
(Xm = 0.5, and 0.9) the opposite happens, namely that the 
smoothed PDF distributions become less Gaussian for any 5Tb 
value. The PDF distributions with and without self-regulation have 
similar shapes, but the presence of low-mass sources makes the 
distribution much less wide. For x m = 0.5 the very brightest 
peaks are fewer than a Gaussian would predict, but there are many 
more intermediate-brightness (5mK < 5Tb — 5Tb < 12 mK with 
no self-regulation, 5Tb — 5Tb < 5 mK with) ones. At the late 
stages of reionization (x m — 0.9) both cases show many more 
bright peaks (10 mK < 5Tb — 5Tb) than a Gaussian would pre- 



dict, although the self-regulated case yields fewer very bright peaks 
(15 mK < 5Tb — 5Tb) than either the corresponding Gaussian or 
the non-self-regulated case. Finally, regardless of the above differ- 
ences in the PDFs, their equivalent widths are very similar for the 
two simulations. 

The PDFs for our fiducial self-regulated high- and low- 
efficiency cases LI, L2, SI and S2 are shown in Figures|26land|27l 
Unlike the presence and self-regulation of low-mass sources pre- 
sented above, which influenced the PDFs significantly, neither the 
source efficiencies nor the box size have any dramatic effect of the 
resulting PDFs. The smaller boxes do not capture well the bright 
wing of the distribution because the highest density peaks are rare 
and the volume in these cases is too small to capture them. The ef- 
fect of varying source efficiency manifests itself by yielding more 
bright peaks during the early stages of reionization and somewhat 
brighter peaks at its middle stages. 

Finally, the results with a varying source models are shown 
in Figures [28] and [29] Here for clarity we just show a representa- 
tive sub-sample of our full simulation suite. Upon inspection sev- 
eral general trends become clear. If only high-mass sources are 
present (model S5) the distributions are noticeably wider, with a 
long non-Gaussian tail at high differential brightness temperatures 
{5T b - 5T b > 30-40 mK) than in the fiducial self-regulation 
case, SI. Conversely, there are many fewer regions with low, but 
positive (i.e. still mostly neutral) differential brightness tempera- 
tures (5Tb — 5Tb < 15 mK). The reason for this is that the mas- 
sive sources form only at the highest density peaks, leaving neu- 
tral many other density peaks which have not yet collapsed. The 
high-density, neutral gas in those peaks is reflected in the non- 
Gaussian tail at high differential brightness temperatures. Lastly, 
model S4 (low efficiency sources and no suppression) yields inter- 
mediate PDF between the fiducial run and the high-mass sources 
only runs. The PDF's are therefore mostly dependent on which 
population of sources is active (high or low mass), but are not very 
sensitive to the details of the reionization history (models S5 and 
S8, not shown here, which have the same source population active, 
but with different efficiencies over time and thus different reioniza- 
tion history yield very similar distributions). These trends are inde- 
pendent of the smoothing employed, as can be seen in Figure [29] 
although naturally the range of differential brightness temperatures 
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Figure 29. Same as in Fig. 1281 but for boxcar smoothing of 5 h 1 Mpc. 




Figure 30. (top) The evolution of the rms of the 21-cm fluctuations, and (bottom) evolution of the skewness of the 21-cm PDFs for (left panels) simulations 
LI (blue), and L3 (black); (middle panels) simulations LI (blue) and SI (green); and (right panels) simulations SI (green), S4 (cyan) and S5 (light red). Shown 
are the results at full simulation resolution (solid lines), smoothed with 3' Gaussian beam and 440 kHz bandwidth (dotted lines) and smoothed with a Gaussian 
beam corresponding to a 2.5km maximum baseline and 440 kHz bandwidth (dashed lines). 



is much reduced by the smoothing. The only new feature found in 
the smoothed data is the non-Gaussian tails for negative 8Tb — 5Tb 
at early times (x m — 0.1). These are result of the H II regions 
(8Tb — STf, < 0) in the non-self-regulated cases growing relatively 
large quickly. Consequently, even at these early times their sizes 
become comparable to the smoothing window size, which results 
in the non-Gaussian tails. These were not present in the cell-wise 
PDF distribution, as the individual cells tend to be either fully ion- 
ized (STt — 8Tb < —25 mK) or mostly neutral. 

The level of non-Gaussianity of the PDF distributions can, 
to a first order, be characterized by their skewness, which in turn 
can be used to d istinguish and extract the reionization signals 
( Marker et al . 2009). In Fig.[30]we show the evolution of the skew- 
ness vs. frequency for selected models. We show the skewness for 
the cell- wise PDF, as well as smoothed with 3' Gaussian beam and 
440 kHz bandwidth and smoothed with Gaussian beam correspond- 
ing to a 2.5km maximum baseline and 440 kHz bandwidth (bottom 
panels). Both sets of beam and bandwidth smoothing are roughly 
as expected for the LOFAR array. We also plotted the rms of the 
correspondingly-smoothed 21-cm differential brightness tempera- 
ture fluctuations (top panels). We note that because of the different 



beam- and bandwidth smoothing employed here the rms values are 
slightly different from the ones shown in Figs. |17|[T9l and|20l 

In all cases, regardless of the specific reionization scenario the 
skewness evolution for the unsmoothed (1-cell) PDFs follows sim- 
ilar pattern. It says at a roughly constant, positive value throughout 
the evolution, until it shoot s up just before o verlap, very similar to 
the behaviour observed in lHarker et al.l 120091) . The the beam- and 
bandwidth-smoothing of the PDFs introduces a significant feature 
in the skewness, whereby it becomes negative during the interme- 
diate stages of the evolution. Interestingly, this dip of the skewness 
to negative values closely corresponds to the rise and peak of the 
differential brightness temperature rms fluctuations, preceding it 
slightly in time. This feature is universal, observed for every reion- 
ization scenario and source model we consider here and suggests 
an interesting approach for a detection and/or independent confir- 
mation of the rise and peak of the 21-cm rms fluctuations during 
cosmic reionization. 

On the other hand, the skewness of the 21-cm PDF distribu- 
tions proves to be fairly insensitive to the source model or reion- 
ization scenario, resulting in only slight changes in the values. The 
skewness of the smoothed PDFs also proves largely independent of 
the box size and resolution and of the details of the interferometer 
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beam assumptions (i.e. if it is fixed in angular size or evolves with 
frequency). This suggests that while the feature in the evolution of 
the skewness is a good indicator of the rise in patchiness, it most 
likely cannot be used for constraining the properties of the ionizing 
sources. 



6 SUMMARY 

In this work we have used a large set of cosmological structure for- 
mation and reionization simulations in an attempt to gain insight 
into what can be learned about the properties of the reionization 
sources based on their observational signatures. In particular, we 
are interested in determining what observations could be used to 
discriminate between certain source models, thereby restricting the 
available parameter space. Here we primarily focused on the red- 
shifted 21-cm signatures, as these can in principle probe the full 
reionization history and offer a wide range of different probes, from 
the mean history, through fluctuations measures like rms evolution 
and power spectra, to PDFs and higher-order statistics which can 
detect non-Gaussian features. 

The observable features of the epoch of reionization derive 
from the gradual mean transition of the IGM from neutral to highly 
ionized state, as well as from the patchiness of that transition. The 
mean transition depends largely on the overall number of ionizing 
photons emitted by the source per unit time, with some correction 
due to recombinations which is position-dependent due to density 
spatial variations. The patchiness, on the other hand depends, in a 
complicated way on the abundances, clustering and efficiencies of 
the ionizing sources. 

Our structure formation simulations confirm previous results 
by us and other groups that the high-redshift halo mass functions 
are inconsistent with either of the widely-used Press-Schechter and 
Sheth-Tormen analytical fits. In particular, the abundance of rare 
halos is strongly underestimated by PS, but over-estimated by ST. 
We find that the nonlinear halo bias is extremely high and very 
scale-dependent. Linear bias regime is only reached at very large 
scales, k > 0.1. For the rarest halos (3-cr and above) linear bias 
regime is never reached even within our largest, 163 Mpc volume. 
Therefore, proper account for the nonlinear bias of halos is impor- 
tant and any calculations assuming linear bias are underestimating 
the halo clustering significantly. 

The Jeans mass filtering of low-mass halos in ionized regions 
and the related self-regulation of reionization results in signifi- 
cantly more extended reionization history and higher integrated 
electron scattering optical depth (by Ar cs ~ 0.01) compared to the 
high-mass source-only scenario with the same overlap redshift, al- 
beit both optical depths are still within the current constraints from 
WMAP Even more significant are the changes in the reionization 
geometry, resulting in corresponding differences in the reioniza- 
tion observables. The 21-cm fluctuations are lower at all scales and 
their PDF distributions are somewhat more Gaussian, although sig- 
nificant non-Gaussianity remains. 

In all our simulations reionization occurs inside-out, with the 
high-density regions being reionized on average earlier than the 
mean and low-density ones. This inside-out nature of reionization 
results in the mass-weighted IGM photoionization rates being con- 
siderably larger, by factor of a few, than volume-weighted ones. 
This should always be taken into account as it can easily skew ob- 
servational measurements depending on the mean density of the 
regions being probed. 

The skewness of the 21-cm PDF distribution smoothed over 



LOFAR-like window shows a clear feature correlated with the rise 
of the rms due to patchiness. This feature does not exist in the un- 
smoothed data, indicating that it is related to the non-Gaussianity 
of the large-scale patchy distribution of 21cm emission. The feature 
exists for any reionization scenario and ionizing source properties 
and thereby provides a different approach for detecting the rise of 
large-scale patchiness and an independent check on other detec- 
tions. 

The peak position of the 21-cm rms fluctuations depends sig- 
nificantly on the beam and bandwidth smoothing size as well as on 
the reionization scenario. As a consequence, it does not always oc- 
cur at 50% ionization fraction as sometimes is claimed, but instead 
can happen for ionized fractions as low as 30% and as high as 70%. 

The simulation volume has only a modest effect on the re- 
sults as long as the typical size of the ionized patches is smaller 
than the volume. However, the fluctuations at large scales (above 
approximately a fifth of the boxsize) are severely affected. This is 
especially important at late times, when the ionized patches grow 
very large. Therefore, at least ~ 100//i Mpc boxes are required to 
model the fluctuations during the late stages of reionization. 

The ionizing source efficiencies and their correlation proper- 
ties introduce clear signatures in the reionization observables. As 
a direct consequence of that, one cannot model low-mass, unre- 
solved sources by simply assigning their emissivity to the resolved 
higher-mass sources as the latter have abundance which is highly 
variable over time and different clustering properties from lower- 
mass sources which provide the bulk of the ionizing photons (see 
Appendix |Al. 

When self-regulation is present there are only minor differ- 
ences between the 21-cm observational signatures resulting from 
high- and low-efficiency ionizing sources, apart from an overall 
shift of the reionization history. The corresponding PDF distribu- 
tions are also very similar, which suggests that the source efficien- 
cies in such models can only be constrained by the overall timing 
of the mean reionization history. 

Scenarios where low-mass sources are completely absent, e.g. 
somehow rendered sterile, are relatively easily distinguishable from 
the ones where they are present (even if strongly suppressed). On 
the other hand, our results suggest that numerous low-efficiency 
sources (case S4) can mimic the effects of suppression (SI). Such 
scenarios therefore might be difficult to distinguish solely based on 
power spectra and similar measurements. However, they might still 
be discriminated through 21-cm PDFs as the no suppression case 
creates many more high-brightness peaks. Similarly, a high-mass 
source only scenario (S5) gives quite similar results to the high- 
efficiency, no suppression case (S3) at the same stages of reioniza- 
tion (albeit these cases overlap at different times for the parameters 
we have chosen), although they do differ in terms of power at large 
scales and in number of bright peaks. 

The results presented in this work should not be considered 
to be an ultimate, realistic prediction of the reionization signals. 
While the assumptions about source efficiencies and suppression 
we made for our fiducial cases are reasonable based on our current 
knowledge and likely bracket the realistic range, the uncertainties 
are still substantial. As more observational data becomes available 
over time it can be used to restrict the parameter space further and 
help us refine our theoretical models, which in turn will provide 
a valuable tool for interpreting the meaning of the observational 
results in terms of early structure formation, source efficiencies, 
suppression mechanisms, etc. In this framework our current study, 
which evaluates in a controlled way the effects of a set of widely 
different assumptions about the sources of ionizing radiation is a 
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very useful step towards a more complete understanding of early 
galaxy formation and feedback. 
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Figure Al. Effective efficiency factors g 7 vs. redshift for simulations S8 
and S9, denned so as to ensure the same total number of ionizing photons 
emitted per atom as in our fiducial simulation, S 1 . 



APPENDIX A: REIONIZATION BY RARE, MASSIVE, 
VARIABLE-LUMINOSITY SOURCES 

Simulations 53Mpc.uvS_le9 and 53Mpc.uvS.lel0 (S8 and S9 in 
Tabled are investigating the effects of keeping the global, volume- 
averaged emissivity of ionizing photons per unit time fixed at each 
redshift, while raising the minimum source mass, assumed to be 
1O 9 M for S8 and 1O 1O M for S9. The overall number of pho- 
tons emitted at each timestep are set to be exactly equal, at all 
times, to the one yielded by our fiducial case, SI. The resulting 
effective efficiencies g 1 are shown in Figure lATI Simulation S8 has 
the same high-mass halo population as our fiducial simulation, but 
no active low-mass sources (LMACHs) at all. Therefore, the ef- 
fective efficiencies start very high, at several thousand photons per 
atom, as the relatively few high-mass sources at early time have to 
'compensate' for the more numerous low-mass sources present in 
the fiducial case which are missing here, as well as for all photons 
emitted before z ~ 19 in our fiducial simulation, during which 
time there are no active sources larger than 10 M@. However, as 
the number of high-mass sources rises exponentially, <? 7 , e // drops 
precipitously, to less than 40 by 2 = 12.6 and less than 20 by 
2 = 11. Towards overlap g 7 , e // settles on ~ 8.7, the value adopted 
in our fiducial case, as by then all low-mass sources are suppressed 
and the high-mass sources are identical to the ones in the fiducial 
case. Note that although in this case the sources belong to the same 
halos as in S5 and overlap is reached at a similar redshift, this S8 
case is different in assuming the same step-by-step total emissivity 
as in our fiducial case, LI, which naturally makes them variable in 
time, unlike case L5 discussed before, which had a fixed photon 
emissivity per unit halo mass. 

Case S9 is still more extreme, since only quite massive halos, 
with masses above 1O 1O M0 are allowed to be active sources. The 
first such massive halos form in our simulation only at 2 = 13.2 
and they remain relatively rare (~ 3 — a) even at overlap (2 = 
8.2). As a consequence, their effective efficiency is very high at 
all times, starting at over 20,000 and reaching ~ 40 at overlap. In 
order to avoid hyper-luminous sources during the first timestep, we 
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Figure A2. Spatial slices of the ionized and neutral gas density from our radiative transfer simulations with boxsize 53 Mpc at box-averaged ionized fraction 
by mass x m ~ 0.50. Shown are the density field (green) overlayed with the ionized fraction (red/orange/yellow) and the cells containing sources (dark/blue). 
Shown are cases S8 and S9. 
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Figure A3. The evolution of the rms fluctuations of the 21-cm background, for beamsize 3' and bandwidth 0.2 MHz and boxcar filter vs. frequency (left) and 
vs. mean mass-weighted ionized fraction (right). Shown are simulations SI (green, short dashed), S8 (dark red, long dash-short dash) and S9 (brown, dotted). 



distributed the photons that were emitted at z > 13.2 in the fiducial 
case over the first several timesteps of run S9. Clearly, both S8 and 
S9 scenarios are not very realistic physically, given this vast range 
of change in the source efficiencies. 

Since in cases S8 and S9 we imposed the same global inte- 
grated ionizing photon emissivities per timestep as in our fiducial 
case SI (but higher minimum source mass), the averaged global 
reionization histories of those two cases closely follow the one of 
the fiducial simulation once the first haloes above the respective 
minimum cutoff form in our volume. The only remaining differ- 
ence is that at early times (z > 11) the ionized fraction in case 
S9 is a little lower than in the other two cases, as a consequence of 
our imposition of a bit more gradual initial release of photons in this 
case, in order to avoid hyper-luminous sources, as explained above. 
However, unlike SI, both simulations S8 and S9 yield x m /x v ~ 1, 
since in those latter cases the ionized patches produced by the few, 



luminous sources present are far less correlated with the underlying 
density field (Fig. |A3t . This is a consequence of the I-fronts quickly 
escaping into the nearby voids, which compensates for the expo- 
nential rise of the number of ionizing sources forming at the high 
density peaks (although we note that even in this case reionization 
remains inside-out, as the ionized regions are still over-dense on 
average). This results in H II region distributions which are clearly 
distinct from the rest. As the cutoff mass increases, there are ex- 
ponentially fewer ionizing sources, which consequently are much 
more efficient (cf. Fig. lAU . Hence, those hyper-luminous sources 
produce correspondingly large H II regions, which are less corre- 
lated with the underlying density field and are more spherical than 
in the other cases, as they are produced by few, but highly clustered 
sources. 

This distinct H II region geometry of cases S8 and S9 also 
yields very characteristic 21-cm signatures. The massive, rare, 
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Figure A4. 21 -cm differential brightness temperature fluctuation power spectra for varying source models. Shown are the epochs at which the ionized fractions 
are (left) x m = 0.1, (middle) x m = 0.5 and (right) x m = 0.9. All cases are labelled by color and line-type, as follows: SI (green, short-dashed), S8 (dark 
red, long dash-short dash) and S9 (brown, dotted). 



highly efficient sources quickly produce very large H II regions and 
thus high rms fluctuations at large scales and a very broad peak, 
with an almost constant value ((ST£) 1/2 = 10.39 - 10.94) for a 
wide range of mean ionized fractions by mass, x m — 0.21 — 0.58. 
There is also no initial dip of the rms fluctuations, which normally 
occurs when the highest density peaks are ionized, but the H II 
regions are still much smaller than the smoothing beam size. In 
models S8 and S9 the ionized patches grow so fast that their typ- 
ical sizes are of order or large than the beam at all times. Such a 
scenario therefore yields a signal which is both stronger and quite 
different from the others. The results for the lower minimum source 
mass cutoff case with same reionization history, S8, show similar 
properties to S9, namely a broad and relatively high rms peak and 
no initial dip. However, the peak value in this case, at ~ 8 mK, is 
noticeably lower than for model S9 (~ 11 mK) and is more similar 
to the typical values for the majority of cases (~ 5 — 8 mK). The 
corresponding PDF distributions (not shown) are noticeably wider, 
with a long non-Gaussian tail at high differential brightness temper- 
atures (8Tb — 5Tb > 30 — 40 mK) than in the fiducial self-regulation 
case, SI. 

The corresponding 21-cm power spectra (Figure lA4l for cases 
S8 and S9 also show significantly higher signal on all scales com- 
pared to our fiducial case S 1 . During the early stages of reioniza- 
tion the power at large scales (k <J 0.4 h/Mpc) for S9 is almost an 
order of magnitude higher than for SI. During the late stages the 
difference decreases considerably, but still remains ~ 2 on large 
(k < 2 h/Mpc), as well as very small (k > 8 h/Mpc) scales. As 
could be expected, the case with lower minimum mass, S8, is inter- 
mediate between S 1 and S9, but much closer to S9 throughout the 
evolution. 

Our scenario S9 is similar to the h i gh mi nimum source mass 
case, S4, considered in iMcOuinn et al.l ( 120071) . These authors set 
the minimum source mass to 4 x 1O 1O M0, somewhat higher than 
in S9. Their ionizing photon production is similarly set to repro- 
duce, step-by-step, the one of their fiducial case. Their results are 
qualitatively similar to what we find. The rare, efficient and strongly 
clustered sources yielded 21-cm power spectra which were higher 
and flatter than in their fid ucial case, with the di fference decreas- 
ing over time (cf. Fig. 17 in lMcOuinn et alj 2007). However, some 
quantitative differences remain, due to the somewhat different ap- 
proach we have taken, as well as some numerical and resolution 



differences. Ap a rt from the higher source mass cutoff adopted by 
IMcOuinn et alj {2007), which results in a stronger source bias, 
other important differences include lower resolution of their N- 
body and radiative transfer simulations, and lack of Jeans mass 
filtering. Unlike our high-resolution simulations, which resolve all 
atomically-cooling halos (M > 10 8 Mp), the N-bod y structure for- 
mation simulations used by IMcOuinn et al.1 ( 120070 resolved only 
halos with mass above 10 9 M©, with lower-mass sources included 
in some cases by sub-grid modelling. More importantly, their fidu- 
cial case (whose photon production per timestep was the basis for 
their high-mass cutoff case S4) yielded late overlap and included 
no Jeans mass filtering (several of their other simulations included 
it, but not this one). Therefore, their photon production per unit 
source mass was necessarily very low, making their fiducial case 
more similar to our low-efficiency case S4 than to our fiducial simu- 
lation S 1 . Finally, we take account of peculiar velocity when calcu- 
lating more precise 21-cm p ower spectra (including redshift space 
distortions) (Mao et al. 201 1 ). Despite these differences, our results 
agree reasonably well on a qualitative level. 



We note that models S8 and S9 are rather unrealistic, as they 
assume unphysically high and time- variable luminosities, as well as 
the suppression of all sources with mass below 10 10 A/q (or, less 
aggressively, 10 9 M© for case S8), for which no clear mechanism 
exists. We have included these models here primarily in order to 
demonstrate, under controlled circumstances, the effect of higher 
source-mass cutoff on the 21-cm observables. Such a higher source 
cutoff mass occurs numerically in simulations with large volumes 
and limited dynamic range (e.g. iBaek et alj |2009j; iThomas et al.1 
2009), and, therefore, it is important to evaluate the level of re- 
liability of such models. Our results show that including only the 
high-mass sources can result in over-estimating the 21-cm rms fluc- 
tuations by up to a factor of 2, while P(k) at small k where the first 
generation of observations will probe, could be over-estimated by 
as much as an order of magnitude at the 50% ionized epoch. It can 
also yield quite a different evolution, even for the same boxsize, nu- 
merical resolution and the same integrated photon emissivity over 
time. One should therefore be aware of these potential pitfalls and 
adjust their modelling accordingly. A better simulation approach 
would be to add the lower-mass, unresolved sources by sub-grid 
modelling. 
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Figure A5. (left) (bottom panel) Mass-weighted reionization histories for cases SI, S6 and S7, each with different Jeans mass filtering threshold, (top panel) 
Ratio of the corresponding mean mass- weighted and volume ionized fractions, x m /x v . (bottom panels) Number of ionizing photons emitted by all active 
sources (thick lines) in the computational volume per timestep and (top panels) cumulative number of photons per total gas atom released into the IGM. 
Vertical lines mark the overlap redshift in each case. All curves on both left and right are labelled by color and line-type, as follows: SI (green, short-dashed), 
S6 (light blue, long dash-short dash) and S7 (red, dotted). 




Figure Bl. Spatial slices of the ionized and neutral gas density from our radiative transfer simulations with boxsize 37 h" 1 Mpc , all at box-averaged ionized 
fraction by mass x m ~ 0.50. Shown are the density field (green) overlayed with the ionized fraction (red/orange/yellow) and the cells containing sources 
(dark/blue). Shown are cases SI, S7, and S6. 



APPENDIX B: DEPENDENCE ON THE JEANS 
SUPPRESSION THRESHOLD 

Here we consider variations of our source suppression threshold, 
with the goal establishing the robustness and validity of our fidu- 
cial model, SI, where we use ^threshold = 0.1. We have ran two 
additional models, S6 and S7, whereby we raise this ionization 
threshold for low-mass source suppression to ^threshold = 0.9 and 
0.5, respectively. In our fiducial case the suppression criterion for 
partially-ionized cells is more aggressive than in these new cases. 
The reasonable value to be adopted for this suppression threshold 
is still quite uncertain at present, thus it is important to check the 
sensitivity of our results to variations in its value. In fact, it is most 
likely that a sharp on-off condition like this is an oversimplification 
and in reality the suppression boundary is gradual, with full sup- 



pression of the smallest galaxies, partial one for intermediate-mass 
galaxies, up to no suppression at all for sufficiently massive galax- 
ies. However, given the current uncertainties, the range of possible 
suppression models is very large and it is difficult to fully explore 
numerically. Instead, we have chosen to consider three very differ- 
ent cases covering the full range of the threshold value, in order to 
evaluate the effect of these uncertainties on the reionization history 
and observables. We note, however, that we consider our original 
source suppression criterion to be well motivated, for the following 
reasons. Although for numerical reasons our suppression criterion 
is defined in terms of ionized fraction, physically it is related to the 
temperature state of the IGM, for which the ionization state is used 
as a proxy. When a given region is photoionized, its temperature 
rises to ~ 10 4 K, with the exact value dependent on the intensity 
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Figure B2. 21-cm differential brightness temperature fluctuation power spectra for varying source models. Shown are the epochs at which the ionized fractions 
are (left) x m = 0.1, (middle) x m = 0.5 and (right) x m = 0.9. All cases are labelled by color and line-type, as follows: SI (green, short-dashed), S6 (light 
blue, long dash-short dash) and S7 (red, dotted). 
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Figure B3. The evolution of the mean 21-cm background (left) and its rms fluctuations for Gaussian beamsize 3' and bandwidth 0.2 MHz and boxcar frequency 
filter (right) vs. observed 21-cm frequency. Shown are simulations (green, short-dashed), S6 (light blue, long dash-short dash) and S7 (red, dotted). 



and spectrum of the ionizing radiation, ranging from ~ 20, 000 K 
for Pop. II stellar spectra and QSO's to - 30, 000 - 40, 000 K for 
Pop. Ill dShapiro et alj|2004l) . This rises the gas pressure and thus 
the Jeans mass, to 1O 9 M0 or more. In order for the low-mass halos 
to be able to re-form in a previously-ionized region its tempera- 
ture should decrease to well below 10 4 K. However, in the mostly 
metal-free gas during these early epochs there is no efficient radia- 
tive coolant available and therefore the main cooling mechanisms 
are the local adiabatic expansion and Compton scattering of CMB 
photons. Since both of these processes are relatively slow and inef- 
ficient we expect that our fiducial more aggressive ionized fraction- 
based source suppression criterion is more physically realistic than 
the milder suppression of the new cases. However, given the sig- 
nificant uncertainties of the Jeans filtering process, which can only 
be properly modelled by hydrodynamical simulations with detailed 



and realistic microphysics, we consider all of these very different 
suppression criteria and study their consequences below. 

The reionization histories and cumulative numbers of ioniz- 
ing photons emitted derived for the three suppression criteria are 
shown in Figure lA5l (left). Clearly, only a very weak suppression 
(^threshold = 0.9, case S7) yields any significant differences. Com- 
pared to our fiducial case SI, many fewer low-mass sources are 
suppressed in S7, and of these a significant fraction are allowed to 
become active again shortly after suppression (since in absence of 
radiation recombinations quickly bring the neutral fraction back up 
above 10%). The evolution of the number of photons produced in 
simulation S6 is up to ~ 2 higher in the middle stages of reioniza- 
tion, while the corresponding number for S7 is essentially the same 
as in S 1 throughout the evolution. 

The Jeans mass filtering nonetheless still has a significant ef- 
fect, keeping the ionized fraction well below the corresponding one 
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for the no-suppression case S3 (cf. Fig. [9] right panels). Eventu- 
ally, the fully-ionized fraction of the volume becomes sufficiently 
large to suppress almost all low-mass sources even with this mild 
suppression criterion and the reionization process slows down until 
sufficient number of massive sources form and are able to finish this 
process and reach overlap. In contrast, the intermediate case, S7 
(^threshold = 0.5) shows only modest differences from the fiducial 
model S 1 , manifesting themselves mostly in bringing reionization 
forward by Az ~ 0.4, compared to Az ~ 1.2 — 3.2, and very 
different shape of the reionization history for case S6. Similarly, 
the integrated electron-scattering optical depth for the mild sup- 
pression case S6 is r cs = 0.111, much higher than in the fiducial 
case (r cs = 0.080), while for the intermediate suppression case 
the increase is much more modest, at r es = 0.089. The cumulative 
number of photons per atom at overlap, ~ 2, is very similar in all 
three cases. 

The variations in the geometry of reionization (Fig. IBlt are 
mostly found in the small-scale structures. There are significantly 
fewer such structures in the weak suppression case S6. The merged 
H II regions are typically slightly larger, as well as rounder and 
with smoother boundaries compared to the fiducial simulation SI. 
Once again the intermediate case S7 is very similar to SI, with 
only minor differences in small-scale features. These visual impres- 
sions are further confirmed by comparing the 21 -cm power spectra 
(Figure lB2l l. The weak suppression case S6 has significantly more 
power at intermediate and large scales (k <C 5) during all stages of 
reionization, more so at late times, but less power on small scales 
than our fiducial case S 1 . On the other hand, the intermediate model 
S7 matches S 1 fairly closely, except for having less power on very 
small scales. 

The 21 -cm mean differential brightness temperature (Fig. |B3l 
left) for the weak suppression case S7 shows an initial steep de- 
cline around v ~ 100 MHz, followed by a sudden change of slope 
at v ~ 125 MHz and a very slow decrease thereafter. This be- 
haviour is quite different from models SI and S6 (which again fol- 
low almost identical evolution), as well as from all other models 
discussed earlier. The only model with a similarly sharp decrease 
of the mean brightness temperature is the high efficiency, no sup- 
pression case S3, which however does not have the same long slow 
evolution tail at late times due to lack of suppression. 

Finally, the differential brightness temperature rms fluctua- 
tions (Fig. |B3l right) for the weak suppression case, S6, peak much 
earlier, at v ~ 110 MHz (but still significantly later than the no sup- 
pression case S3, which underlines the importance of even a very 
weak low-mass source suppression) reaching (ST 2 ) 1 ' 2 ~ 8 mK. 
Uniquely, this model exhibits a very long tail of slow decline of 
the rms fluctuations beyond the peak. This is related to its very dif- 
ferent (and, as we argued earlier, possibly less physically realistic) 
suppression model. Case S6 also exhibits significant fluctuations 
in the differential brightness temperature fluctuations, once again 
indicating that this suppression model might be less physically re- 
alistic than our standard suppression model. The intermediate sup- 
pression model S6 follows the same evolution as the fiducial case 
SI, but shifted to slightly earlier time. 

In summary, all results prove fairly insensitive to the precise 
value of the Jeans suppression threshold assumed, as long as it is 
not at the very weak suppression limit. Both ^threshold = 0.1 and 
0.5 yield essentially the same evolution, apart from a slight off- 
set in time. On the other hand, a very high suppression threshold 
(^threshold = 0.9, i.e. weak suppression) results in a very different 
(and somewhat unstable) evolution with several characteristic ob- 
servational signatures. We have, however, argued above that such a 



weak suppression is likely less realistic physically. Our suppression 
model therefore proves quite robust to a threshold variation within 
the plausible range. 
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